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TECHNICAL PUBLICATION 


SHOCK STRUCTURE ANALYSIS AND AERODYNAMICS IN A WEAKLY 

IONIZED GAS FLOW 

1. INTRODUCTION 


Unlike the fully ionized case, the structure of a shock wave propagating through a weakly ionized 
gas is not well understood. 1 Interest in the weakly ionized case, however, has been rekindled by experi- 
mental research suggesting that the performance of a supersonic missile or aircraft can be significantly 
improved by weakly ionizing the air through which the associated shock wave propagates. 2-8 

Since the ion/neutral scattering cross section is large and the electron mass is negligible, one may 
assume that ions and neutral particles in the weakly ionized gas have the same temperature and velocity 
and that the electron momentum is negligible. Based on these assumptions together with a quasi-neutral- 
ity assumption, Appleton and Bray developed a “heavy particle/electron” model for a weakly ionized gas 
with both elastic and inelastic (ionizing) collisions, 9 while McIntyre, et al. used a derivative of this model 
to characterize the relaxation and recombination processes in a weakly ionized gas. 10 

Under the same assumptions, Magretova, et al. and Liu and Glass developed heavy particle/elec- 
tron models for an ionizing shock wave propagating through a non-equilibrium weakly ionized gas. 11 ’ 12 
Simulations using the Magretova, et al. model at high temperatures in the Mach 26-29 range, show a small 
electron temperature peak and a jump in both the charged and neutral particle densities at the shock, and 
simulations using the Liu and Glass model at room temperature and Mach 13 yield similar results. 11 ’ 12 
(See ref. 1 for a review of the Magretova, et al. model and simulations. 1 ’ 11 ) Popovic and Vuskovic and 
Saeks, et al. dropped the quasi-neutrality assumption in their heavy particle/electron models (though still 
requiring that the electron velocity equal the heavy particle velocity). 13 ’ 14 In addition to the effects indi- 
cated above, they saw a large electric field, on the order of 10 5 -10 7 V nr ! , at the shock wave and a reduc- 
tion in the shock intensity, which is consistent with the experimental observations. 2-8 

The object of this Technical Publication (TP) is to analyze the structure of a shock wave propa- 
gating through a weakly ionized gas and characterize its influence on the aerodynamics and performance 
of a missile or aircraft. Lor this purpose, an electrofluid dynamics model composed of six classical con- 
servation laws and Gauss’ Law was developed, while the remaining conservation laws are implied by 
the assumptions that u l = u n>{) and T l = T n and that the electron momentum is negligible. Viscosity is 
included in the model to correctly define the shock thickness, but quasi-neutrality is not assumed since the 
spatial scale of the shock structure may be less than a Debye length. As such, both the electric field and the 
electron drift velocity are included in the model. 


This model is divided into a fluid dynamics model and an electrodynamics model. The fluid 
dynamics model is derived in section 2. 1 using the multispecies fluid flow model developed by Sutton and 
Sherman with the addition of the viscosity terms. 15 Since the electron velocity may be several orders of 
magnitude greater than the heavy particle velocity, an analysis of this model reveals that the electron body 
force term in the energy equation (equivalently the electron ohmic heating term) is the primary mechanism 
through which the electrodynamics interact with the fluid flow. As such, the role of the electrodynamics 
model is to provide sufficient electrodynamics equations to evaluate the body force/ohmic heating term in 
the energy equation. The resultant electrodynamics model, formulated in section 2.2, is made up of three 
individual species conservation laws plus Gauss’ Law and the required auxiliary equations. The source 
and sink terms used in the electrodynamics model are summarized in section 2.3. 

Although the electrofluid dynamics model is formulated in vector notation, the primary goal of 
this TP is to analyze the structure of a normal shock wave propagating through a weakly ionized gas and 
identify the underlying physical processes. As such, it suffices to solve a one-dimensional version of the 
electrofluid dynamics model, restricted to the region a few tens of micrometers upstream and downstream 
of the shock where the ionization interacts with the fluid flow. Indeed, in this region, the fluid dynamics 
model reduces to a pair of initial value problems upstream and downstream of the shock derived in sec- 
tion 3.1, which are connected by analytic intrashock and shock width models derived in section 3.2. The 
flow chart for a one-dimensional shock structure analysis algorithm based on the resultant model appears 
in section 3.3. 

A detailed analysis of the structure of a Mach 4 shock wave at an altitude of 40,000 ft using a 1-eV 
source is presented in section 4. 1 together with a discussion of the physics underlying the key features of 
the shock structure. Additional shock stmcture simulations appear in appendix A. The normally intuitive 
concepts of flow and shock intensity are formalized and used to summarize the properties of a series of 
shock structure simulations — as a function of Mach number, altitude, and the source temperature— in sec- 
tion 4.2. Furthermore, it is shown that the shock structure effectively bifurcates into two cases: (1) Where 
the ionization does not significantly effect the shock structure and (2) where the ionization reduces the 
shock to its minimum possible intensity, with a small transition region in between. 

An analysis of the flow behind a shock wave propagating through a weakly ionized gas is devel- 
oped in section 5. Given the bifurcation of the shock structure, the two cases are considered separately 
using a classical analysis for the case where the ionization does not significantly effect the shock structure, 
and a modified analysis for the case where the ionization reduces the shock to its minimum possible inten- 
sity. The required modified pressure and density ratio formulae are developed in section 5.1 and used to 
formulate a set of oblique shock curves for the case where the ionization reduces the shock to its minimum 
possible intensity in section 5.2. The modified oblique shock curves are used to determine if the shock 
is attached or detached and to compute the shock angle when it is attached, and an analysis of the flow 
behind the shock for the attached and detached cases is performed in section 5.3. 

In section 6 the results of an aerodynamic analysis of a two-dimensional hypersonic lifting vehicle 
in which ionization is used to reduce the shock intensity to its minimum possible value, comparing its 
performance at various Mach numbers and angles of attack to an identical vehicle without ionization, are 
considered. 
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Finally, in section 7, the conclusion, the assumptions underlying the electrofluid dynamics model 
are reviewed, the shock structure analysis and simulations are compared with experiment, and the possible 
applications of the technology in aeronautics are commented upon. 
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2. ELECTROFLUID DYNAMICS MODEL 


To formulate the electrofluid dynamics model, consider a fluid composed of neutral gas mol- 
ecules, single positively charged ions, and electrons where each particle species is characterized by the 
variables summarized in table 1. Here the superscript, s, indicates the particular species ( s = n, i, e for 
neutral molecules, ions, and electrons, respectively), whereas no superscript is used for the entire ion- 
ized gas treated as a single fluid. Furthermore, dual superscripts are used as required: s = ed to denote the 
electron drift velocity, s = ie to denote variables that apply to both the ions and electrons, etc. Clearly, 
p s = m s n s and M s = u s / C\ while it follows from the ideal gas law that p s =n s kT s , where/: is Boltzmann’s 
constant, and e s = p s j[y s _ 1 j p s for a perfect gas where y v is the specific heat ratio for the given par- 
ticle species. 15 ’ 16 In the following analysis a weakly ionized gas; i.e., n e , n l <sc n 11 , is considered with no 
externally imposed electric or magnetic fields. An electric field due to charge separation may, however, be 
present over a distance of a few Debye lengths. 17 For nonrelativistic motion, however, the force due to a 
magnetic field, |±w s / C pjeH , is much less than that due to an electric field, eE, and as such, any magnetic 
fields due to charge transport may be ignored. Furthermore, u l = u n and T l = T n is assumed since the ions 
and neutrals have approximately the same mass and collision cross section and, therefore, flow together. 
See reference 17 for a detailed analysis of these effects, while a summary of all assumptions used in the 
derivation of the fluid dynamics model is provided in table 2. 


Table 1. Variables and constants characterizing each particle species in a weakly ionized gas. 


Variable 

Symbol 

Units 

Variable 

Symbol 

Units 

Velocity vector 

u s 

m s -1 

Number density (concentration) 

n s 

ITT 3 

Pressure 

P s 

N m -2 

Internal energy per unit mass 

e s 

J kg- 1 

Density 

P s 

kg m -3 

Mass 

m s 

kg 

Temperature 

F 

K 

Speed of sound 

C s 

m s _1 

Mach number 

M s 

1 

Body force 

f s 

N rrr 3 

Entropy 

S s 

J kg- 1 K" 1 

Mass flow 

<P 

kg m -2 s _1 

Electric field 

E 

V m _1 

Current density 

f 

Am -2 

Electric potential 

<D 

V 

Ohmic (Joule) heating 

0 s 

J rrr 3 s _1 

Shock width 

A 

m 

Mean free path 

l 

m 

Dynamic viscosity 

'] 

kg rrr 1 s _1 

Second viscosity 

$ 

kg m _1 s _1 

Heat addition 

Q 

J nr 3 s _1 

- 

- 

- 

Constant 

Symbol 

Units 

Constant 

Symbol 

Units 

Electric charge 

e 

C 

Boltzmann's Constant 

k 

JK” 1 

Permittivity of free space 

e 0 

C 2 N- 1 nr 2 

Speed of light 

c L 

m s _1 

Mass of a neutral molecule 

m n 

kg 

Mass of an electron 

m e 

kg 

Ionization potential 

£ i 

J 

Ionization temperature 

01 

K 

Specific heat ratio 

yS 

1 

- 

- 

- 
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Table 2. Assumptions and bounds used in deriving the electrofluid dynamics model. 


Physical Assumptions 

Approximation 

Implied Equality 

Justification 

The neutrals and ions have the same velocity 
and temperature. 

n , i 
U =U 

ft rjn 1 

The neutrals and ions have (approximately) the same mass and 
collision cross section and therefore flow together. See ref. 1 7 for a 
detailed analysis of this effect. 

There are no externally imposed electric or 
magnetic fields 


The model is intended to analyze the effect of a weakly ionized gas 
on a shock wave in a free space environment without externally 
imposed fields. 

Internally generated magnetic fields are 
ignored. 

H = 0 

For non-relativistic motion the magnetic force, (±u s IC L )eH , is much 
less than the electric force, eE, allowing H to be ignored. 17 

Steady state analysis 

o 

II 

T’O' 

The model is intended for steady state applications. 

Perfect gas 

p s =n s kT s 

The model is intended for aerodynamic applications in air (or other 
perfect gases). 

Ions are singly ionized 

Z'= 1 

Since the gas is weakly ionized, the probability of multiply ionizing a 
single molecule is small. 

The are no radiative heat losses 

2 = 0 

Thermal radiation from the shock is negligible because of its small 
width. 18 

Parameter Bounds 

n e < lO^ 4 Xn" n' < 10 -3 Xn” T e <l0 2 xT n u e <10 5 xu n 


Unlike the ions, the electrons will be accelerated by an electric field, where present, with their velocity 
limited by collisions with the heavy particles. Indeed, u e may be expressed in the form u e = u n + u ed + u ei 
where u ed is the electron drift velocity due to an electric field and u er is the random electron velocity due 
to scattering. In general, the electron drift velocity, u ed , can be modeled by a differential equation charac- 
terizing the acceleration of the electrons by the electric field less the momentum lost in electron neutral 
collisions. 17 In this application, however, the approximate expression, u ed = -eE 1 / m e v en , suffices, where 
v en is the electron/neutral collision frequency. Typically, u ed is two or more orders of magnitude greater 
than u n , whereas the random electron velocity may be one to two orders of magnitude greater than u ed . 1 7 
Since the random electron velocity is due to the scattering of the electrons off of the neutrals, its direction 
is random, and as such, the mean random velocity of the electrons is zero, allowing the modeling of the 
mean electron velocity by 


eE 


e en 
m v 


(i) 


Note that in the presence of a positive flow, u n = u > 0, u e must also be positive. Indeed, if 
u e (x) < 0, at some point, x, the electrons would flow upstream while the ions would flow downstream. The 
resultant increase in upstream electron density would, however, cause E(x) to decrease via Gauss’ Law; 
i.e., become more negative, which would increase u e to a positive value. As such, a negative electron 
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velocity in the presence of a positive flow is unstable. Moreover, it follows from the above expression for 
u e , that E is bounded above by E < u n m e v en /e, which is approximated by E < 0 since even a small positive 
E will produce a “weak electron trap”. 

The electrofluid dynamics model is divided into two components, (1) A fluid dynamics model 
developed in section 2.1, characterizing the fluid flow properties of the entire weakly ionized gas treated 
as a single fluid; and a electrodynamics model developed in section 2.2. The source and sink terms used 
in the electrodynamics model are summarized in section 2.3. 

2.1 Fluid Dynamics Model 

Following Sutton and Sherman, 15 define the fluid dynamics variables for the entire weakly ionized 
gas treated as a single fluid by letting 


P=p n +P'+P e 


and 


n , i , e 
p = p + p + p 


be the sum of the individual species densities and partial pressures, respectively, and 


u — 


p n u n +p'u l + p e u e 


( 2 ) 


(3) 


(4) 


and 


£ . pv yv+py (5) 

P 

be the mass density weighted average of the individual species velocities and internal energies, respec- 
tively. Given the relatively small weight on the electrons in the fluid dynamics variables, they typically 
track the corresponding neutral particle variables within 1 percent, given the bounds for a weakly ionized 
gas indicated in table 2. This is not to imply that the ionization has no effect on the neutral particle flow 
(sec. 4), but rather that its effect on the neutral particle flow tracks its effect on the entire fluid. As such, 
one can use the fluid dynamics variables for the weakly ionized gas to approximate the neutral particle 
variables in the electrodynamics model developed in section 2.2. 

Using the above defined fluid flow variables and following Sutton and Sherman with the addition 
of the viscous terms, 15 use the following set of fluid flow equations: 
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V • (pu) = 0 


(6) 


p{u • V)m + Vp - 7?V 2 w - (C + 77 / 3)' V(V • m) = /" +/ + / e 


(7) 


and 


V* 


f 

f 

i 2 A A 



u 


p 


— — + £ 

u 



2 


V 
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/ 


+ V* pu=f n •u n + f' •u l + f e • u e +Q 


( 8 ) 


to model the weakly ionized gas. Here, 77 and £ are the dynamic and second viscosity coefficients, 
respectively; f s ,s= n, i, e is the body force on the three particle species and Q is the heat added to the 
ionized gas. Note that, since the above equations define a fluid dynamics model for the entire weakly ion- 
ized gas treated as a single fluid, the interspecies collision processes are implicit in the equations, and do 
not require explicit models. 


In the momentum and energy equations (7) and ( 8 ) we have expressed the body force,/*, in terms 
of the individual particle species rather than the entire weakly ionized gas, since the electric field, which is 
the source of the body force in the present problem, effects each particle species differently. Letting/* be 
the body force per unit volume due to the electric field on each of the three particle species in the energy 
equation ( 8 ), f n *u n = 0 , while 


f e »u e = - en e E *u e = -en e u e • E = j e »E = £2 e (9) 

reduces to the ohmic (or Joule) heating of the ionized gas due to the body force on the electrons, where 
j e = - en e u e is the electron current density, and similarly, f l »u = j l »E -Q l . Since the electron drift 
velocity is two or more orders of magnitude greater than the ion velocity, in the presence of an electric 
field, the ion ohmic heating can be neglected in comparison to the electron ohmic heating. As such, the 
body force terms in the energy equation reduce to f *E = £l e . Note that the ohmic heating term, £2 e , only 
represents that part of the energy imparted to the electrons by the electric field, which is lost to the neu- 
tral particles via collisions, with the remainder going into the random electron velocity (see ref. 17 for a 
detailed analysis of the ohmic heating process). 

Unlike the body force terms in the energy equation, the body force terms in the momentum equa- 
tion are not multiplied by velocity, and as such, their effect on the fluid flow variables is normally neg- 
ligible compared to £2 e . Moreover, since the electric field only exists over a distance of a few Debye 
lengths on either side of the shock if one integrates in the direction normal to the shock over the interval 
\x up , v c | n J where E is non-zero : 17 
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\ XA \f l {q) + f e {q) dq= \ X An l (q)-n e (q) eE(q)dq 

Am-* l— — A l— — 

U P up 


1 


dn 


up 

£ 0 z?2 


dE , , 
£ o— (?) 

dx 


Km 

dE 2 (g) 

up 

dx 


dq 


f z\q) 


dn 

'up 


: ^ [^(^dn) £(*up)] 


= y[0-0] = 0 , 


( 10 ) 


where the equality \n l - n e je = e^dEjdx is the differential form of Gauss’ Law in the direction normal 
to the shock. As such, the net first order effect of the ion and electron body force terms in the momentum 
equation is zero — sometimes called the washout effect— while f n = 0 . Although higher order and/or 
local effects are not precluded, equation (10), together with the relatively small magnitude of the body 
forces compared to Q e , is sufficient justification to neglect the body force terms in the momentum equa- 
tion.* 


Finally, letting e = ( p n e n + p'e 1 + p e e u j p~ p/(y — l)p and substituting equations (9) and (10) 
into equations (6)— (8) yields the required fluid dynamics model for the ionized gas treated as a single fluid, 
as follows: 

• Continuity equation: 


V»(pw) = 0 . 


( 11 ) 


• Navier Stokes equation: 


p(w*V)« + Vp-77V 2 M-(£-i-77/3)V(V»u) = 0 . (12) 

• Ionized gas energy equation: 


V* 


( ( 

P 


\ \ 


- + - 


(y-i)p 


+ V • pu=Q t 


(13) 


* Since the electric field produced by charge separation across the shock only extends a few Debye lengths, 17 the number of electron / neutral collisions in the 
interval [x up , x dn ] may be insufficient to transfer the momentum imparted to the electrons by the electric field to the entire weakly ionized gas. In this case, 
the electron body force term drops out of the momentum equation (eq. (7)), negating the “washout effect.” Since fscO, however, the remaining ion body 
force term in the momentum equation is additive to the ohmic heating term in the energy equation (eq. (8)), enhancing the interaction between the ionization 
and the fluid flow variables. As such, the “washout analysis” yields a conservative approximation in this case. 
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Although the heat addition term in the energy equation (eq. (8)) can be used to model the effects 
of an external heating source on the shock, our goal in the present paper is to characterize the effects of a 
weakly ionized gas on the shock, independently of thermal effects, while following Lowke and Murphy, 
thermal radiation from the shock is negligible. 18 As such, set Q = 0 in equation (13). 


2.2 Electrodynamics Model 

Given the fluid dynamics model of equations (1 1)— (13), it remains to add the appropriate electro- 
dynamics equations required to evaluate Q e . To this end, define an electrodynamics model composed of 
the following equations. 

• Electron velocity equation: 


e eE 

u = u 

e en 
m v 


Electron heating equation: 


-n e k(u e • Vjr e + n e kT e (v»u e ^ = 77 + 0-12'’ + <7 


Ion continuity equation: 


V • [nu - D 


'Y(V)) = n' e 


Electron continuity equation: 


V.|nV-D e V(n e )| = n ie 


Poisson’s equation: 


V 2 <£ = — —in 1 - n e ) 

e 0 V ’ 


(14) 


(15) 


(16) 


(17) 


(18) 


and 


E = -V0 . 


(19) 
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Electron ohmic heating: 


Q e =-en e E »u e , (20) 

where the approximation u l ~u n ~u is used to eliminate the species variables, u n and u' from the elec- 
trodynamics model. 

The electron heating equation (eq. (15)) is the energy equation for the electrons without the kinetic 
energy term. 15 Here, Q e is the electron ohmic heating term defined in equation (9), I] is the thermal 

energy lost in elastic collisions, using a standard “hard sphere” model, and 0 is the thermal energy lost 
in inelastic; i.e., ionizing, collisions. 15 Finally, <7 denotes the ionization source in the electron heating 
equation. For the present purposes, rather than attempting to model a realistic source, choose <7 to main- 
tain a prescribed minimum electron temperature in the region around the shock where the ionization 
interacts with the fluid flow, T e (x)> T ^ rc , jc up < x < x dn , where T^ rc is typically on the order of 1 eV. 
Detailed models for t? and 0 are developed by Sutton and Sherman, Hoffert and Fien, and summarized in 
section 2.3. 15 ’ 19 


In the charged particle continuity equations (eqs. (16) and (17)), n e - n ie+ -n ie ~ characterizes 
the ionization/recombination rate. For the ionization rate, n e+ , a single step ionization model is used, 
where ionization is due to an electron colliding with a neutral and ionizing it from the ground state. 15 ’ 17 
For the recombination rate, h ie ~ , a model developed by Stevefelt et al., which combines radiational and 
collisional recombination rate terms with an empirical correction term, is used. 20 Finally, D 1 and D e are 
the ion and electron diffusion parameters, respectively, though D l may be replaced by the ambipolar ion 
diffusion parameter. 17 These models are developed by Sutton and Sherman and Raizer and summarized 
in section 2.3. 15 ’ 17 


The electric field, E, is expressed as the gradient of a scalar potential, <£>, which satisfies Poisson’s 
equation, in equations (18) and (19), while the electron velocity, and electron ohmic heating equations 
(eqs. (14) and (20)) were derived in section 2.1. 

Each of the terms used in the electrodynamics model, as well as the collision frequencies, v en and 
v m , and mean free path that appear throughout the model can be expressed in terms of the fluid dynamics 
variables, p, p, and u and/or the individual species parameters in the electrodynamics model. As such, the 
fluid dynamics and electrodynamics models of equations (1 1)— (20) constitute a self-consistent electro fluid 
dynamics model. Ignoring the explicit substitutions of equations (14) and (20), the resultant model is 
composed of one vector and six scalar equations in one vector and six scalar ( u , p, p, T e , n l , n e , and &) 
unknowns. 
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2.3 Source and Sink Models 


The ionization/recombination rate term, n ie = h ie+ - h ie ~ , in the charged particle continuity equa- 
tions (eqs. (16) and (17)) may expressed as the difference of the ionization rate, n ie+ , and the recombina- 
tion rate, n e . For the ionization rate, use the single step ionization model, 


■ie+ — ^ n 
n — — Ci 

3 


£j€ 


r £ilkT e 



8xlC )~ A kT e 


Km 


W 2 , 





( 21 ) 


where ionization is due to an electron colliding with a neutral and ionizing it from the ground state. This 
model is derived from the models in references 15 and 17, where £/ is the ionization potential of a neutral 
molecule and Cj is an empirical constant, 14 ’ 17 while the n l term in equation (21) is added to enforce the 
single ionization hypothesis. For the recombination rate model, use 


■je- 


1.55xl0“ 16 (r e ) 


- 0.63 


+ 3.61x10 ~ ll \T e 


— 2.18 / \ 0.37 

K 


+3.8xl0“ 21 (r e ) 


- 4.5 


e i 
n n 


( 22 ) 


developed by Stevefelt, Boulmer, and Delpech, which combines radiational and collisional recombination 
rate terms with an empirical correction term. 20 

In the electron heating equation (eq. (15)), I] is the thermal energy lost in elastic collisions, using 
a standard “billiard ball” model, 15 


77= 3 — v en k(T-T e ) , (23) 

m n 1 ' 

where v en is the electron/neutral collision frequency. 17 Similarly, (j) = - ri ie+ £j is the thermal energy lost 
in inelastic; i.e., ionizing, collisions. For an electron to ionize a neutral molecule from the ground state, 
£j Joules are transferred from the electron to the molecule. Note that the ionizing energy is radiated 
by the molecule when it recombines and relaxes to the ground state and is normally not returned to the 
electrons. 

Finally, D l and D e in the charged particle continuity equations are the ion and electron diffusion 
parameters 


D 


kT 


m n v ul 


and D e 


kT e 


e en 
m v 


(24) 
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where v is the ion/neutral collision frequency, though D l may be replaced by the ambipolar ion diffu- 
sion constant D a = kT e / m ll v m P 
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3. SOLUTION OF THE FLUID DYNAMICS EQUATIONS ACROSS A SHOCK WAVE 


The purpose of this section is to develop a procedure for analyzing the structure of a shock wave 
propagating through a weakly ionized gas using the electrofluid dynamics model developed in section 2. 
Of course, one could simply solve equations ( 1 1)— (20) in one dimension, but this would require a priori 
knowledge of the downstream boundary conditions and a full downstream solution of the fluid dynamics 
equations to capture the shock. Alternatively, one can analyze the shock structure without explicit knowl- 
edge of the downstream boundary conditions by assuming that the shock wave exists and (without loss 
of generality) that it is located at x us = 0 while ignoring the effect, if any, of the downstream boundary 
conditions on the intrashock flow. Indeed, this approach is implicit in the classical case where one uses 
the Rankine-Hugoniot equations to compute the pressure and density ratios across the shock, combined 
with analytic shock width and intrashock flow models, if required, to analyze the shock structure without 
explicit reference to the downstream boundary conditions. 16 


Over the past few years, several authors have modified the Rankine-Hugoniot equations to include 
the effect of ionization on the pressure and density ratios across the shock wave. 21,22 In practice, however, 
the effect of the ionization on the shock structure is spread out over a region a few tens of microns upstream 
and downstream of the shock as indicated in figure 1. As such, Murray and Saeks proposed to replace the 
single control volume used to define the modified Rankine-Hugoniot equations with a sequence of control 
volumes upstream and downstream of the shock, starting with the given upstream flow and chaining the 
pressure and density ratios across each control volume to compute the pressure and density ratios across 


the entire region where the ionization interacts with the flow. 21 

Upstream Side of Shock, x us — . 

\ 


. — Downstream Side of Shock, x ds 

/ 

Upstream Flow 

1 


Downstream Flow 

1 



Electric Field 


Figure 1. One dimensional geometry for shock structure analysis. 


In the present analysis, this process is taken to the limit, letting the thickness of the control vol- 
umes go to zero with the chain of pressure and density ratios converging to initial value problems in the 
regions upstream and downstream of the shock connected by an analytic intrashock model. The initial 
value problems in the regions upstream and downstream of the shock, connected by an inviscid (zero 
thickness) shock are formulated in section 3.1, while the inviscid shock is replaced by a viscous intra 
shock model and an associated shock width model in section 3.2. Finally, a computational algorithm for 
analyzing the structure of a shock wave propagating through a weakly ionized gas, without reference to 
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the downstream boundary conditions, is obtained by combining these techniques with a one-dimensional 
version of the electrodynamics model in section 3.3. 


Although the electrofluid dynamics model using vector notation is defined in section 2, for the 
actual shock structure analysis the one-dimensional geometry used for this TP is illustrated in figure 1. 
The flow is moving forward along the x-axis, which is translated so that the upstream side of the shock is 
located at x us . The subscript “us” is also used to denote the fluid flow and electrodynamic variables at the 
upstream side of the shock; p us = p(x u s ) , E us , etc. Similarly, the downstream side of the shock is located 
at x ds , with the subscript “ds” used to denote the variables at the downstream side of the shock. Unlike 
a classical shock, the ionization can affect the flow upstream of the shock via electron and ion diffusion 
across the shock. As such, let x up be a point upstream of the region where the ionization begins to interact 
with the flow, and use the subscript “up” to denote x and the variables at x up , and similarly for a point x dn 
and the variables at x dn , downstream of the region where the ionization ceases to interact with the flow. As 
such, the structure of a shock wave propagating through a weakly ionized gas may be fully characterized 
by the fluid flow and electrodynamic variables in the interval [x up , x dn ] . 

3.1 Modified Rankine-Hugoniot Analysis 

To formulate the initial value problems that subsume the Modified Rankine-Hugoniot equations in 
the regions upstream and downstream of the shock where the ionization interacts with the fluid flow, start 
with a one-dimensional version of the fluid dynamics equations (eqs. (1 1)— (13)): 

— (p«) = 0 , (25) 

dx 
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du dp , d^u 

pu — + Tj ' — y 

dx dx dx 


0 


(26) 


and 


d_ 

dx 


( ( 



V V 


(y-i)p 


\ a 

u 

2 J 


+ — {pu) = Q e 
dx y ’ 


(27) 


where 77 ' = (£ + 4 77 / 3) . To simplify the analysis, treat the viscosity term in equation (26) as an external 
force rather than an intrinsic part of the momentum equation by computing the second derivative using 
data from a previous iteration of the shock structure analysis algorithm or smoothing and extrapolating 
previously computed data points on the present iteration (sec. 3.3). In either case, this process decouples 
the second derivative from the unknown velocity, allowing us to treat it as an independent variable, which 
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2 2 

is denoted by d u / dx . Although this process will smooth out some higher order effects, it is sufficient 
for the present purpose. 


dp p d 

From equation (25), it follows that the mass flow, (p = pu is a constant, and — = u . Using 

dx u dx 


the equality cp-pu to eliminate p as a variable in equations (26) and (27), yields 


dp 

dx 



+ n 



and 


d 


(..2 


u pu 


2 (y-l)<p 


+ ^-(pu) = Q e , 
dx 


while expanding equation (29) yields 


du p du u dp du dp _ e 


(pu — + 


- + - 


+ p — + u £2 , 


dx (y-l)Jx (y-1 ) dx dx dx 


(28) 


(29) 


(30) 


and 


du YP du yu dp _ ^ e 


(pu 1- 


- + - 


dx (y- 1) dx (y- 1) dx 


Substituting equation (28) into equation (31) yields the series of equalities 


du YP du yu 
(pu — + — u —. + ■ 
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£2 e 
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(32) 


(33) 
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(r- 1 ) 


(yp-(pu) C ^ t - = n e 
dx 


2 

yu ,d u 
— - — n 

(7-1) dx 2 ’ 


(34) 
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and the desired expression for 


du _ (r~l) r, e Y“ 

dx ( yp-(pu ) ( yp-(pu ) dx 2 


while substituting equation (35) into equation (28) yields 
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dx 
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Q e 
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,d 2 u 

dx 2 


(35) 


(36) 


Equations (35) and (36) define a set of “fluid state” equations that can be used to compute p(x) 
and u(x), given Q e , d 2 u t dx 2 , and initial values for p and u, and one can express the remaining fluid 
flow variables in terms of p and u. In practice, however, the “fluid state” variables M and C are prefered; 
i.e., the Mach number and speed of sound of the flow. Indeed, with the aid of the equalities u = CM and 
p = (pC / yM , one can transform equations (25) and (26) into an equivalent set of “fluid state” equations, 

tfC_ (r~l) , y(y-l) M 3 A (37) 

dx 2 (p C 2 (p dx 2 

and 

<im _ (r-i)( l+ r M2 ) M | y (2+(7 -i)m 2 )m% , g 

dx lip |i-A/ 2 ) C 2 2<p |l-A/ 2 ) C dr 2 

in M and C. Finally, use Sutherland’s equation to compute 77 , 23 Indeed, in this application temperatures 
less than 200 K are rarely encountered, and as such, Sutherland’s formula may be approximated by 
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(39) 
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Furthermore, in this application, the “fluid state” equations are used in the regions upstream and 
downstream of the shock but not in the shock. As such, it is not expected to encounter extreme compres- 
sion or expansion waves, allowing the approximation of the second viscosity by the dynamic viscosity. 16 
Setting rj' = £ + ri/3~lri/3, yields 


7 a \m n 

77' *-x 1.458x10 “ 6 xJ— xC = aC , (40) 

_J7 

where a = 1.69 x 10 . Substituting equation (40) into equations (37) and (38) and setting u = CM now 

yields the desired form for the “fluid state” equations, 


dC & + aY{Y- 1) M 3 C d 2 CM 

dx 2cp \-M 2 C 2(p [\-M 2 \ dx 2 


(41) 


and 


dM _ (y-1) 1 

[l + yM 2 ] 

* M Q e ay! 

[2 + (y-l)M 2 

) M ^ d 2 CM 

dx 2 (p 

1 -M 2 

C 2 2 (p 

1 

[l -M 2 ) 

dx 2 


which can be solved for C and M given Q e \ d 2 u / dx 2 and initial values for C and M. Moreover, given the 
“fluid state” variables, C(x) and M(x), one may compute the remaining fluid flow variables via 

w = CM; p = (p!CM\ and p = (pC/yM . (43) 

The most significant feature of equations (41) and (42) is the singularity at M- 1, which is a mani- 
festation of the physical fact that the flow cannot pass through Mach 1 continuously. 

Indeed, since Q e > 0 , the Mach number of the upstream supersonic flow is reduced by the ohmic 
heating process, eventually producing a shock wave where it jumps across the singularity at Mach 1 to a 

9 

subsonic value. Since the sign of ( 1 - M ) changes at Mach 1, the Mach number of the resultant subsonic 
flow then begins to increase under the influence of the ohmic heating process. Unlike the supersonic flow, 
however, the subsonic flow cannot jump across the singularity, and is accelerated toward Mach 1 by the 
singularity in equation (42), with the rate of acceleration damped by the viscosity term. 16 Moreover, if the 
flow reaches Mach 1, the sign of ( 1 - M ) changes, trapping it at Mach 1. This results in a Mach number 
curve typical of that shown in figure 2 (a), while the pressure curve follows an inverse trajectory illustrated 
in figure 2 (b). In general, whenever the flow is supersonic, M and u decrease in the ionization, while p and 
p increase, and vice-versa when the flow is subsonic. See section 4. 1 for a detailed analysis of a Mach 4 
flow in 40,000-ft air that has been weakly ionized by a 1-eV generator. 
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(a) x 10- 5 (b) x10“ 5 

Figure 2. (a) Mach number and (b) pressure across a Mach 4 shock at 40,000 ft with a 1-eV source 


3.1.1 Solution of the Fluid Dynamics Equations Across a Shock Wave 

The solution of the fluid dynamics equations across a shock wave is as follows: 

(1) Given Q e and cl 2 u / dx 2 , initialize equations (41) and (42) at v up with the upstream flow 
variables, M up , C up , and solve (41) and (42) for M{x ) and C(x) in the interval [x up , -TisJ 5 where x up 
is a point upstream of the region where the ionization interacts with the flow, and x us is the point at the 
upstream side of the shock. 

(2) Use the (classical) Rankine-Hugoniot equations to compute the flow variables at a point, 
x t |s = x us + 8 , at the downstream side of an inviscid shock, given the fluid flow variables at x us , where 8 
is a single step in the numerical process used to compute the flow variables. 16 

(3) Initialize equations (41) and (42) with M ds and C ds at x ds , and solve (41) and (42) for Mix) 
and C(x) in the interval [x ds , .x dn ] , where .t dn is a point downstream of the region where the ionization 
interacts with the fluid flow. 

This process yields a complete solution of the fluid dynamics equations in the region around the 
shock where the ionization interacts with the flow, excepting the inviscid shock, which will be replaced by 
a viscous intrashock model in section 3.2. 
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3.2 Intrashock Model 


Although an inviscid fluid dynamics model producing an ideal (zero-thickness) shock is used in 
classical shock structure analysis, 16 ’ 24 in this application a viscous intrashock model is required to prop- 
erly model the electron energy losses in the shock. To this end, the inviscid shock used in section 3.1 is 
replaced with an analytic intrashock model derived from the Navier-Stokes equation. 16 

Following Landau and Lifshitz, in the one-dimensional case, 16 if the inviscid fluid dynamics model 
indicates the presence of a shock in the interval between v us and v c | s , with pressure p us at the upstream 
side of the shock and p ( [ s at the downstream side of the shock, one approximates the viscous effects in the 
shock by interpolating p mv in the interval [v us , v t | s ] with a hyperbolic tangent, 




(Pds-Ais) tanh 


7.X X us A'gjg 

2A 


+ {Pds Pus ) 


’ ^us < x ^ Lis 


p vts (x) = p inv (x) ; *<* us , x ds <x , 


(44) 


and similarly for p vis (v) and u vls (v) . Here, the shock width, A, is typically a few mean free paths for 

the upstream neutral gas. 16 ’ 24 Note, in the following, the superscripts “mv” and “vis” have normally been 
dropped, distinguishing between the inviscid and viscous models by context. 

Even though A is termed the “shock width” in the literature, 16 ’ 24 it is, in actuality, the slope of the 
hyperbolic tangent at the center of the shock. Indeed, less than half (46 percent) of the pressure (density, 
velocity, etc.) jump across the shock takes place over a distance of A, as illustrated in figure 3, while a 
distance of 5 A is required for 99 percent of the jump. As such, 5 A is used when it is desired to model the 
distance required for “most” of the jump across the shock. 



Figure 3. Effective shock width required to model 99 percent of the jump across the shock. 
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Although Shapiro provides an analytic formula for the shock width, A, it is quite complex. 24 Alter- 
natively, Landau and Lifshitz provide an estimate of A, in the form (ref. 16) 


8a 

P 2 (/Lis-Pus)^ 2 (l/p)/3p 2 )^ 


(45) 


2 / 2 2 \ 

together with estimates for a ~ A / C and (d (l / p) / dp ~ 1 / pp , where A is the upstream mean free 
path for the neutral gas. Substituting these expressions into equation (23) and evaluating p at the center of 
the shock, p = ( p ds + p us ) / 2, yields 

A 4 (p ds + Pus) 

* ~ r(Pds-Pus) ’ ; 


which may be used to compute the shock width. Indeed, if one divides equation (46) by p us , and substi- 
tutes the classical Rankine-Hugoniot pressure jump formula into the resultant expression, 16 equation (46) 
reduces to 


PcL +1 2 7 M u 2 s -(y-l) + i 

^ 4 (Pds +Pus) _ 4 p us _4 (7 + l) 4 yM us + 1 ( 47 ) 

* r(Pds-Pus) 1 7 2yM 2 s -(y-1) 1 y 2 M 2 S -1 

Pus (y+ 1 ) 

yielding an alternative expression for the shock width in terms of the Mach number of the flow at the 
upstream side of the shock. 

The shock widths predicted by equation (47) range between 2 and 3 mean free paths for Mach 
numbers greater than 2 in a monotonic gas, increasing to infinity below Mach 2, and the shock width for 
a diatomic gas is =20 percent greater than the shock width for a diatomic gas. As such, the shock widths 
predicted by equation (47) are consistent with those predicted by Shapiro and Shapiro and Kline up to 
Mach 10, 24 ’ 25 while properly accounting for the difference between monotonic and diatomic gases. 24 
Above Mach 10, however, equation (47) converges to 4 / y, while Shapiro’s model diverges. 

Unfortunately, equation (47) is not applicable in the presence of ionization, since M us does not 
account for the effect of the ionization on the shock width. As such, in the presence of ionization, revert 
to equation (46) to compute the shock width directly from the pressure distribution, p. Indeed, with the 
upstream side of the shock at x us and the downstream side of the shock at v ds = x us + 5 A (using the “99 
percent” approximation of fig. 3), p us = p(x us ) and p ds = p(x us + 5 A), while equation (46) reduces to 

A_4 (p(*us + 5A ) + P(*us)) 

* r(p(*us+ 5 A )-p(*us)) 
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which can be solved for A, given p{x ) . Indeed, since Pds > Pus , the expression on the right-hand side of 
equation (48) is a bounded positive continuous function of A, and the expression on the left-hand side of 
equation (48) goes from 0 to 00 as A goes from 0 to As such, the two expressions must cross for some 
A > 0, guaranteeing that equation (48) has a positive real solution. Moreover, since the viscous and invis- 
cid variables coincide at the upstream and downstream sides of the shock, one can use either p vls or p un 
in equation (48). As such, to incorporate the viscous intrashock model into the solution of equations (41) 
and (42) add the following steps to the solution described in section 3.1.1: 

(4) Solve equation (48) for A using p ,nv . 

(5) Interpolate the fluid flow variables between x us and x ds = x us + 5 A with hyperbolic 
tangents. 


3.3 Computational Process 

Although the electrofluid dynamics model of equations (1 1)— (20) is potentially applicable to a 
wide variety of electrofluid dynamics problems, the development of a complete computational electrofluid 
dynamics code to implement the model is well beyond the scope of the present investigation, where it 
suffices to solve the electrofluid dynamics equations for a one-dimensional flow a few tens of microns on 
either side of the shock. A flow chart for the shock structure analysis algorithm developed for this purpose 
appears in figure 4. In the flow chart, the subscripts “k” and “/’ denote the inner and outer loop iteration 
indices, respectively. The fluid flow variables,/^-, p-, and Uj and the electric field, Ej from the / th outer loop 
iteration are fed back to initialize the (/+l)th outer loop iteration while the electron density from the (j , k)th 
inner loop iteration is fed back to initialize the (j, k+\ )th inner loop iteration. Finally, the computational 
process is initialized with the upstream fluid dynamics variables; p up , p up , and w up ; and “initial values” 
for the fluid dynamics variables, p 0 , p 0 , and w 0 , and the electric field, E {) . Also note that the ohmic heating 
term is the only interface from the electrodynamics model to the fluid dynamics model, and the fluid flow 
variables connect the fluid dynamics model to the electrodynamics model. 
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Figure 4. Shock structure analysis algorithm flow chart. 
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The shock structure analysis algorithm illustrated in figure 4 is composed of the following five 
processes: 

(1) Electron velocity and heating equations. 

(2) Ion continuity equation. 

(3) Electron continuity equation and Poisson’s equation. 

(4) One-dimensional fluid dynamics problem with ohmic heating over an interval a few tens to 
hundreds of microns on either side of the shock wave. 

(5) Analytic intrashock model into the solution of the one-dimensional fluid dynamics problem. 
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4. SHOCK STRUCTURE IN A WEAKUY IONIZED GAS 


The purpose of this section is to summarize the results of a series of shock structure simulations 
using the electrofluid dynamics model developed in sections 2 and 3 together with the iterative algorithm 
summarized in section 3. Section 4.1 provides a detailed analysis of the structure of a Mach 4 flow through 
a normal shock wave in 40,000-ft air that has been weakly ionized with a 1-eV source— additional analy- 
ses appear in appendix A. In section 4.2, the concept of “shock intensity” as the geometric mean of the 
pressure and density jumps across the shock is formalized, allowing the summarization of the results of a 
shock structure simulation in a single plot, and comparison of the results of a large number of simulations 
parameterized by Mach number, altitude, and source temperature. 

4.1 Shock Structure Analysis 

The purpose of this section is to summarize the results of a shock structure analysis for a Mach 
4 flow at an altitude of 40,000 ft with a 1-eV source, while providing some insight into the underlying 
physical phenomena. In the electrofluid dynamics model developed in section 2, the source maintains a 
minimum electron temperature (typically between 0.5 and 1.5 eV) across the region where the ioniza- 
tion interacts with the flow (between x up and x dn in the geometry defined in figure 1). Such a source will 
produce electron and ion fractions on the order of 10 _7 -10 -5 . Since the ions and electrons are created and 
annihilated simultaneously, however, the generator by itself will not produce an electric field and the asso- 
ciated ohmic heating required for the electrodynamics model to interact with the fluid dynamics model. 
Rather, the electric field is due to a jump in the ionization rate at the shock wave, with the additional elec- 
trons diffusing farther upstream than the ions, thereby producing a negative electric field across the shock 
wave and the associated ohmic heating. 

In previous analyses the jump in ionization rate across the shock has been attributed to: 

• A jump in the electron temperature across the shock driven by the differential term in the electron 

heating equation (eq. ( 1 5)): 13,14,26—28 j n inviscid case, where the loss terms in the electron heating 

equation are negligible in the shock, one can show that the electron temperature ratio across the shock 
is equal to the pressure ratio raised to the 2/3 power, producing an orders -of- magnitude increase in the 
ionization rate at the shock. When one incorporates the true shock thickness into the model, however, 
the loss terms in the electron heating equation (primarily the transfer of electron energy to the neutrals 
via ohmic heating) damp out much of the electron temperature jump. 

• The jump in ionization across the shock: 11 ’ 12 Since the ionization in this model is due to electron neu- 
tral collisions, the ionization rate is proportional to the product of the gas density and the electron den- 
sity. 15 ’ 17 As such, when the gas density jumps across the shock the ionization rate increases, producing 
additional electrons that further increase the ionization rate (in a positive feedback mode damped by the 
loss of electrons due to transport and recombination), resulting in a jump in the ionization rate across 
the shock that exceeds the jump in neutral density across the shock. 
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Indeed, in our 40,000 ft Mach 4 example with a 1-eV source, the electron temperature jump across 
the shock is completely damped out by the ohmic heating losses as illustrated in figure 5 (a), while the 
ionization rate illustrated in figure 5 (b) tracks n x n e (figs. 6 (a) and 8 (a)), indicating that it is due to the 
jump in gas density across the shock. In particular, the ionization rate continues to increase, tracking the 
increase in electron density, downstream of the point where the gas density saturates. 
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Figure 5. (a) Electron temperature and (b) ionization rate across a Mach 4 shock 
at 40,000 ft with a 1-eV source 



Although the electron diffusion constant is typically several orders of magnitude greater than the 
ion diffusion constant, 17 when one incorporates the effect of the electric field, with the attractive forces 
between the ions and electrons increasing the ion diffusion length; i.e., ambipolar diffusion, 17 and the 
electron drift velocity decreasing the electron diffusion length, the electron diffusion length is typically 
only a few tens of microns greater than the ion diffusion length, yielding an electron density which is only 
slightly more diffuse than the ion density as indicated in figure 6 (a). This is, however, sufficient to gener- 
ate a negative electric field on the order of 10 5 -10 7 V/m at the shock, as indicated in figure 6 (b). Finally, 
given the electron density and the electric field one can evaluate the electron velocity and the ohmic heat- 
ing term which serves as the interface between the electrodynamics model and the fluid dynamics model, 
as indicated in figures 7 (a) and 7 (b). Here, the sharp drop in electron velocity across the shock is due to 
the increase in the electron/neutral collision frequency, which roughly tracks the increase in pressure, 17 
across the shock. 
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Figure 6. (a) Electron and ion densities and (b) electric field across a Mach 4 shock 
at 40,000 ft with a 1-eV source. 
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Figure 7. (a) Electron velocity and (b) ohmic heating across a Mach 4 shock at 40,000 ft 
with a 1-eV source. 


The primary fluid dynamics variables for our Mach 4, 40,000 ft flow with a 1-eV source are illus- 
trated in figures 2 and 8. 
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(a) x 10 -5 (b) x 10 -5 




(c) x 10 -5 (d) x 10- 5 

Figure 8. (a) Density, (b) velocity, (c) temperature, and (d) speed of sound across a Mach 4 shock 
at 40,000 ft with a 1-eV source. 


Unlike the classical case where the upstream flow is constant, the diffusion of electrons (and ions) 
across the shock (fig. 6 (a)) extends the ohmic heating term upstream of the shock (fig. 7 (b)), where it 
effects the upstream flow. In particular, the pressure and density increase while the Mach number and flow 
velocity decrease upstream of the shock. Fortunately, the decrease in pressure ratio across the shock, due 
to the decreased upstream Mach number, more than compensates for the pressure increase at the upstream 
side of the shock and yields a net decrease in the pressure at the downstream side of the shock. A similar 
analysis is applicable to the density at low Mach numbers. At high Mach numbers, however, where the 
density ratio approaches its limiting value, the drop in density ratio may be insufficient to compensate for 
the increase in density at the upstream side of the shock, yielding a net increase density at the downstream 
side of the shock. 
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In the shock, the ionization further reduces the pressure and density, while increasing the Mach 
number and velocity. This effect was analyzed in detail using the modified Rankine-Hugoniot equa- 
tions. 14 ’ 21 ’ 22 Moreover, the ionization in the shock increases the shock width (eq. (46)). 

The ionization continues to decrease the pressure and density, and increase the Mach number and 
velocity downstream of the shock. This effect may dissipate at some subsonic Mach number (fig. 10 (c)), 
or as in the present example, may drive the Mach number back to 1, where it is “trapped” by the singularity 
in the “fluid state” equations (eqs. (41) and (42)) since a subsonic flow cannot jump across Mach 1. Indeed, 
given the nature of this singularity, as the Mach number approaches 1, the flow is accelerated at an increas- 
ing rate towards Mach 1 . As such, the downstream Mach number is driven back to Mach 1 in many of our 
simulations, while the pressure and density are driven to corresponding minimum values (sec. 5). This 
results in the “overshoot” effect observed in figures 2 and 8, where the fluid flow variables jump across the 
shock to an intermediate value (between their “no ionization” value and their final value), and then drop 
to their final value downstream of the shock. 


The temperature and speed of sound across the shock are plotted in figures 8 (c) and 8 (d), respec- 
tively. Although the ionization reduces both the downstream pressure and density, it reduces the down- 
stream density more than the downstream pressure, thereby causing the downstream temperature to increase 
with ionization and similarly for the downstream speed of sound. Finally, the “kink” in the speed of sound 
curve can be explained by examining the |l- yM 2 jy^l- M 2 j term in equation (41), which changes 
sign at both M = 1 and M = 1 / Jy . As such, if the Mach number drops below 1 / <Jy downstream of the 
shock before returning to Mach 1 (as in this example— see figure 2 (a)), dC/dx changes sign four times. 
Specifically, C increases upstream of the shock where M > 1 , decreases immediately downstream of the 
shock where 1 / Jy < M < 1 , increases while M < \f Jy , and decreases again as M approaches Mach 1 
downstream— 1 / yfy < M < 1 . 


For a point of comparison, selected plots from two additional shock structure analyses are pre- 
sented in figures 9 and 10. The first example in figure 9 is the same as above with the source increased 
from 1 to 1.4 eV, chosen to illustrate the sensitivity of the shock structure to the source temperature. As 
above, the electron temperature spike is completely damped, but the ionization rate increases rapidly with 
the increase in gas density across the shock as illustrated in figure 9 (a). This, in turn, causes the electron 
density, the electric field, and the electron velocity to increase, combining to increase the ohmic heating, 
shown in figure 9 (b), by almost two orders of magnitude. As such, the ionization upstream of the shock 
is sufficient to reduce the Mach number of the flow at the upstream side of the shock to Mach 2, while 
increasing the pressure at the upstream side of the shock to 6xl0 4 Pa, illustrated in figures 9 (c) and 9 (d), 
respectively. This, in turn, increases the Mach number and decreases the pressure in the shock. Finally, 
the ionization in the shock is sufficient to drive the downstream Mach number back to 1 and the down- 
stream pressure to its maximum value of 1.8xl0 5 Pa within the shock, thereby eliminating the overshoot 
observed in the previous example. 
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Figure 9. Shock structure analysis for a Mach 4 shock at 40,000 ft with a 1.4-eV source: 

(a) Ionization rate, (b) ohmic heating, (c) Mach number, and (d) pressure. 


The second example shown in figure 10 is a Mach 12 flow at an altitude of 10,000 ft with a 1-eV 
source, chosen to illustrate some of the effects which occur at high Mach numbers. In this example, the 
downstream pressure jump due to the high Mach number and the associated increase in the electron/neu- 
tral collision frequency, 17 reduces the downstream electron velocity by an order of magnitude, as illus- 
trated in figure 10 (a). This, in turn, reduces the downstream ohmic heating, shown in figure 10 (b), to a 
point where it is insufficient to drive the pressure to its minimum value, as indicated in figure 10 (c), and 
similarly for the density (not shown). Moreover, in this case, the electron temperature losses due to ohmic 
heating are insufficient to fully damp out the jump in electron temperature at the shock, as indicated in 
figure 10 (d). 
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Figure 10. Shock structure analysis for a Mach 12 shock at 10,000 ft with a 1-eV Source: 

(a) Electron velocity, (b) ohmic heating, (c) pressure, and (d) electron temperature. 


4.2 Shock Intensity 

The term “shock intensity” is commonly used as a qualitative measure of the pressure and/or den- 
sity jump across a shock wave. The goal of this section is to formali z e this concept, with the goal of using 
it to quantify the effect of ionization on the shock structure without resorting to the full set of plots used 
in section 4.1. To this end, the intensity of a flow at jc is defined as the geometric mean of the pressure 
and density at x , l(x) = p(x)p(x ) . Since p and p increase and decrease together, while u = <p/ p,l(x) 
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simultaneously characterizes the combined effect of the shock on p, p, and u. Moreover, with the aid of 
some algebra, it follows that 


l(x) = jp{x)p(x) 


vVry 


1 

M(x) ‘ 


(49) 


As such, i(x) defines a single parameter that simultaneously characterizes the effect of both the ionization 
and the shock on the primary flow variables. 

Finally, define the intensity of a shock wave as the ratio of the downstream flow intensity to the 
upstream flow intensity: 
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Since the downstream flow is subsonic, the inequality, M dn < 1 , lower bounds i sllock downstream 
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(51) 


implying that there is a limit to the extent to which ionization can reduce the intensity of a shock wave (by 
62 percent at high Mach numbers in air); i.e., if the upstream Mach number of the flow is 3, then the “inten- 
sity ratio” across the shock cannot be reduced below 3 by ionization. Although this appears to contradict 
several experiments where ionization appeared to fully disperse a shock wave, 7 ’ 14 ’ 29 ’ 30 this is believed to 
be an artifact of the Schlieren and shadowgraph diagnostics used in these experiments (sec. 7). 


The normalized flow intensity, i(x)j ^/\ ip p up (which equals the shock intensity at x = x dn ), for 
the 40,000-ft, Mach 4, 1-eV example of section 4.1 is plotted in figure 11, together with the downstream 
lower bound implied by equation (51) and the normalized flow intensity without ionization. 
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Figure 11. Normalized flow intensity for a Mach 4 flow at 40,000 ft with its downstream 

lower bound (lower dashed line) and the “no ionization” case (upper dashed line). 


The normalized flow intensity for a series of simulations at a fixed altitude (40,000 ft), parameter- 
ized by the Mach number and the source temperature appear in figure 12, while the normalized flow inten- 
sity for a series of simulations at a fixed source temperature (1 eV), parameterized by the Mach number 
and the altitude appear in figure 13. Although the downstream flow intensity is M up ; i.e., the minimum 
possible value, in most of the simulations in figures 12 and 13, one can quantify the effect of ionization 
in these simulations by the distance downstream of the shock required for the flow intensity to reach M up . 
The trends observed in these plots and the underlying physical processes are summarized below. 

4.2.1 Source Temperature 

Since the ionization rate can be expressed in the form r[r e ^n e n , where r [r e \ is an increasing 
function of electron temperature, 15 ’ 17 the ionization rate increases when T e is increased at constant Mach 
number and altitude. This, causes n e to increase, further increasing the ionization rate and the electron 
density in a positive feedback mode damped by the loss of electrons due to transport and recombination. 
The resultant increase in electron density causes an increase in the electric field, which in turn, increases 
the electron drift velocity. As such, all three variables that define the ohmic heating process increase, which 
causes Q e to increase cubically and the shock intensity to decrease via the “fluid state” equations (41) 
and (42). As such, the shock intensity decreases with increasing electron temperature. This is consistent 
with the simulations presented in figure 12, while the cubic increase in Q e accounts for the high sensitiv- 
ity of the shock intensity to small increases in electron temperature. 
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Figure 12. Normalized flow intensity versus Mach number and source temperature at 40,000 ft. 
(Note the change in scale on the Mach 8 and 12, 0.8-eV plots.) 
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Figure 13. Normalized flow intensity versus Mach number and altitude with a 1-eV source. 
(Note the change in scale on the Mach 12, 10,000-ft plot.) 
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4.2.2 Mach Number 


When one increases the Mach number at constant electron temperature and altitude, the down- 
stream pressure and density increase. The increase in downstream density causes the ionization rate, 
r{T e )n e n , to increase, producing additional electrons that further increase the ionization rate and electron 
density as indicated above. On the other hand, the increase in downstream pressure increases the electron/ 
neutral collision frequency, which reduces the electron drift velocity and the ohmic heating term. Finally, 
the increase in upstream velocity required to increase the Mach number at constant altitude, increases the 
mass flow, (p , which scales down the ohmic heating term, Q e , in the “fluid state” equations (eqs. (41) and 
(42)), reducing the effect of the ionization on the shock wave. 

Although the downstream density can theoretically increase by a factor of 6 (in air), the ioniza- 
tion reduces the increase in downstream density Indeed for the 40,000-ft 1-eV simulations of figure 12, 
the density ratio across the shock only increases from 1.28 at Mach 1.5 to 1.81 at Mach 12. As such, the 
effect of the increased downstream density can be neglected compared to the effect of the increased down- 
stream pressure and mass flow, which increase from 1.76 to 85 and 136 to 1,085, respectively. As such, 
the increased downstream pressure and mass flow dominate the increase in downstream density, causing 
the shock intensity to decrease with increasing Mach number, which is consistent with the simulations 
presented in figures 12 and 13. 

4.2.3 Altitude 

When one increases the altitude at constant electron temperature and Mach number, the gas density 
and pressure decrease, while the mass flow decreases with the decreasing density. Replicating the above 
analyses, the decreased density decreases Q e , while the decreased pressure decreases the electron/neutral 
collision frequency, increasing the electron drift velocity and Q e . and the decreased mass flow scales up 
the ohmic heating term, Q e , in the “fluid state” equations (eqs. (41) and (42)). Given these contradic- 
tory effects, one cannot predict the effect of an altitude change on the shock intensity, as indicated by the 
simulations in figure 13. 

4.2.4 Summary Conclusion 

Finally, note that the decreased effect of ionization at high Mach numbers (and possibly altitudes) 
does not preclude its application since one can readily compensate for these effects with a small increase 
in the electron temperature given the high sensitivity of ohmic heating process to changes in T e . 

Although the continuity of the electrofluid dynamics model with respect to the underlying param- 
eters is manifested by the distance downstream of the shock required for the flow intensity to reach its 
minimum value in the simulations of figures 12 and 13, the effect of the ionization on the downstream flow 
(x > ) is “almost” binary. 

Indeed, the effect of ionization is either negligible with the downstream intensity roughly equal 
to the “no ionization” case, or the ionization level is sufficient to drive the downstream flow intensity 
to M u p with a small transition region between the two extremes (indicated by the 10,000 ft, Mach 12, 1-eV 


35 


simulation in fig. 13). This is due to the singularity in equations (41) and (42), which literally draws the 
flow intensity to M up , if the Mach number begins to approach to 1. Indeed, our simulations indicate that 
a change in electron temperature of 0.01 eV is usually sufficient to cause the downstream flow intensity 
to jump from one extreme to the other, allowing us to treat the downstream flow as “effectively binary.” 
In the following we refer to the two extreme cases as “with ionization” and “without ionization,” respec- 
tively, ignoring the transitory case. 

Indeed, this effect suggests the possibility of analyzing the impact of ionization on the underlying 
missile or aircraft, without recourse to a full shock structure analysis, by reducing the analysis to the two 
cases, “with ionization” and “without ionization” (sec. 5). Also note that if ionization is used to reduce the 
flow intensity for a normal or detached shock wave, the underlying missile or aircraft must be designed to 
operate in the transonic flow behind the shock wave. 
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5. TWO-DIMENSIONAL AERODYNAMICS 


The purpose of this section is to analyze the air flow through and behind a shock wave propagat- 
ing through a weakly ionized gas. Given the bifurcation of the shock structure indicated in section 4, the 
two cases are analyzed separately, using a classical analysis “without ionization,” and a modified analysis 
“with ionization” developed in this section. Modified pressure and density ratio formulae with ionization 
are developed in section 5.1 and used to formulate modified oblique shock curves in section 5.2. The 
modified oblique shock curves are used to determine if the shock is attached or detached and to compute 
the shock angle when it is attached, “with ionization.” Analysis of the flow behind the shock for both the 
attached and detached cases is performed in section 5.3. 

5.1 Modified Pressure and Density Ratios 


The purpose of this section is to develop an analytical model for the pressure and density ratios 
across the shock with ionization. Begin by defining “pressure and density effective Mach numbers” to 
quantify the effect of ionization by matching the pressure and density ratios with ionization to the pressure 
and density ratios without ionization at a lower Mach number. Recall that the pressure and density ratios 
without ionization are given by the normal shock relations: 16 ’ 31 

Pd„_ 2 yMu 2 P -(r-i) and p d „_ (r+i)M u 2 p (52) 

p »p (r +1 ) Pup (r+i)M„ 2 p +2 ' 


As such, one may define “pressure and density effective Mach numbers” by replacing the down- 
stream pressure and density, p^ n and p ( j n , in equation (52) with the downstream pressure and density, 
/i t [ n and p dn , with ionization, and the actual upstream Mach number, M up , with pressure and density 
effective Mach numbers, M up and M up , satisfying 

£dn = 2r^up-(y-l) md Pdn = (r + l)'V?p t (53) 

Pup (r + 1) Pup (r-i)M u 2 p + 2 


Here an “upper bar” is used to distinguish the parameters with and without ionization. Since 
Pdn < /kin an( l Pdn < Pdn > while the expressions in equation (53) are monotonically increasing, 
1 < M up < A/ up and 1 < M up < M up . Unfortunately, M up and M up do not, in general, coincide (except 
when there is no ionization or M up = 1 ). As such, following the model used to define shock intensity, define 
the effective Mach number of a flow with ionization to be the geometric mean of M up and M up : 

M\x p = JWup^up • (54) 
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Clearly, when there is no ionization, M up = ^M up M up = M up , while based on simulations using the 
structure analysis algorithm with ionization, 


M 


up 


-^up-^up ^ ’ 


(55) 


Although we are not aware of any physical justification for equation (55), numerous simulations were 
performed that ranged from Mach 1.5 to Mach 12, with electron temperatures between 0.6 and 1.4 eV, at 
altitudes between 10,000 and 70,000 ft in which equation (55) is satisfied with an average of 4 percent 
error and a standard deviation of 0.5 percent. As such, equation (55) is used as an empirical model in the 
remainder of the analysis. See appendix B for an approximate but fully analytic derivation of the same 
result. 

Given the bifurcation of the downstream flow indicated in section 4.2, the above equations com- 
pletely define the effective Mach number (except for the transitory case). In particular, M up = ^M up 
with ionization and M up = M up without ionization. 

By combining equation (55) with equations (50) and (53), one may solve for M up and M up in 
terms of M up , with ionization. Indeed, in this case, 
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Now, upon substituting equation (55) (in the equivalent form M up 
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which can be solved for M up and M up , in terms of M up , yielding the expressions 



( 58 ) 
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In combination with the equality M up = M up = M up without ionization, equations (58) and (59) fully 
characterize the effective Mach numbers (except for the transitory case). 


Finally, upon substituting equations (58) and (59) into equation (53), one may express the pressure 
and density ratios across the shock with ionization directly in terms of the actual upstream Mach number 
via 


and . (60) 
Pup (r + 1) Pup yM up + l 

By comparing these expressions with the classical normal shock relations of equation (52), one 
may verify that the pressure and density ratios across a normal shock wave are always reduced with ion- 
ization. Comparing the corresponding temperature ratio 

Pdn Y^up + 1 

?dn _ ^up _ ( y + 1) 

r up £dH (7+l)^up 

Pup ymI v + \ 

with the classical temperature ratio; 16 however, one concludes that the downstream temperature ratio 
across a normal shock wave is increased with ionization, which is consistent with the simulations. 

Finally, one can compare the pressure and density ratios with ionization computed numerically 
using the above described shock structure analysis simulations with the pressure and density ratios pre- 
dicted by equation (60). Indeed, in a total of 43 simulations ranging from Mach 1.5 to Mach 12, with 
electron temperatures between 0.6 and 1.4 eY, at altitudes between 10,000 and 70,000 ft, the normal 
shock relations of equation (60) match the simulated data with an accuracy of 5 percent and 3 percent for 
the pressure and density ratios with ionization, respectively, and standard deviations of 0.7 percent and 
0.3 percent, respectively. 

Although the analytical model of equation (60) was “calibrated” using a single data point from 
the numerical data (eq. (55)), which is indirectly an average of the product of the computed pressure and 
density ratios over the entire array of simulations), its ability to predict the individual pressure and density 
ratios with the indicated accuracy supports the voracity of both models. 
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Finally, note that the possibility of deriving the expression M up = ^M up for the effective Mach 
number from first principles has not been totally exhausted and remains an open question. (Also see 
app. B for an approximate but fully analytic derivation of this expression.) 

5.2 Oblique Shock Curves With Ionization 

Given the bifurcation of the shock structure identified in section 4.2, one can use the classical 
oblique shock curves to determine if a shock wave is attached or detached and to compute the shock angle 
without ionization. 31 The purpose of this section is to develop a set of modified oblique shock curves 
with ionization. To derive the oblique shock curves with ionization, analyze the geometry shown in fig- 
ure 14, where % is the turning (or nose cone) angle of the flow and (3 is the shock angle with ionization. 
Invoking the standard trigonometric analysis used to express the density ratio across the shock in terms of 
X and [3 , 31 yields 


^- = cot(/?)tan(/I- x) ■ (62) 

Ptln 



Figure 14. Geometry for the oblique shock curves. 
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Then equation (58) is used to express the Mach number of the normal air flow through the shock 
wave as a function of the density ratio across the shock with ionization: 
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Finally, express the normal upstream Mach number in terms of the actual upstream Mach number, 
Ad up = M” p / |sin (7j)) , obtaining the desired relationship, 
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(64) 


between the upstream Mach number, the shock angle with ionization, and the turning angle. The resultant 
oblique shock curves with ionization are plotted in figure 15 along with the classical obliques shock curves 
“without ionization.” Note that for any given Mach number the shock curve “with ionization” is always to 
the right of the corresponding shock curves without ionization, demonstrating the broadening of the shock 
wave with ionization, which has been widely observed experimentally. 3 ’ 6 ’ 29 ’ 30 ’ 32 Carried to the limit, this 
implies that in many cases of interest where the shock would normally be attached, the shock with ioniza- 
tion is detached; i.e., the horizontal line for the specified turning angle in figure 15 is above the specified 
Mach number curve. 



Figure 15. Oblique shock curves with ionization (solid) 
and without ionization (dashed). 
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5.3 Downstream Flow 


This purpose of the section is to analyze the flow downstream of a shock wave with ionization. To 
this end, the modified oblique shock curves developed in section 5.2 are used to determine if the shock 
is attached or detached and to compute the shock angle when it is attached. The flow through an attached 
shock and around a two-dimensional wedge is analyzed in sections 5.3.1 and 5.3.2. Unlike the classical 
case without ionization, where the shock is attached for most realistic turning angles and Mach numbers, 
with ionization the shock is often detached for medium and large turning angles and/or low Mach num- 
bers. As such, the flow through a detached shock and around a two-dimensional wedge is analyzed in 
section 5.3.3. 

5.3.1 Attached Shocks With Ionization 

The analysis of the flow through an attached shock and around a two-dimensional wedge with 
ionization is similar to the classical analysis of the flow through an attached shock and around a two- 
dimensional wedge with the exception of the pressure and density ratios across the shock, 31 and the fact 
that the Mach number of the downstream flow normal to the shock is always 1 (eq. (51)). The resultant 
flow pattern is illustrated in figure 16, where an “upper bar” is used to denote variables that change with 
ionization. As in the classical case the attached shock is also wedge shaped with a shock angle [3 deter- 
mined from the oblique shock curves with ionization (fig. 15), or by solving equation (54) for [3 , given 
the wedge angle % an< J the upstream Mach number A/ up . Moreover, it follows from equation (64) that a 
streamline crossing the shock at a point x, will turn by an angle of (3 - % as it passes through the shock, 
and then flow parallel to the wedge. 31 


\ 



Figure 16. Flow through an attached shock and around a two-dimensional wedge with ionization. 
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To compute the pressure and density, p & n and p dn , in the streamline at the downstream side of 
the shock, one invokes the modified pressure and density ratios with ionization of equation (60), using the 
Mach number of the upstream flow normal to the shock, M” p = M up sin(/j) , in lieu of the actual upstream 
Mach number, obtaining 

Pin rKp si "(i )) 2+1 md p in (r + l)KpSin(i))~ (65) 

Pu P ( 7 + 1 ) Pu P y(M up sin(^))“ + 1 


Furthermore, since the angle between the shock and the streamline downstream of the shock 
is P - % , while the Mach number of the downstream flow normal to the shock with ionization is 1, 
M dn = l/ sin(/i - x) • 

Finally, since the flow in the streamlines downstream of the shock and parallel to the wedge is 
neither compressive or expansive, it follows from Bernoulli’s equation that the above flow parameters, 
evaluated immediately behind the shock, are constant over the entire downstream portion of the stream- 
line. 16 Moreover, since the shock angle is independent of the point where the streamline crosses the shock, 
the downstream flow parameters are same for all streamlines. As such, the downstream flow is uniform 
over the entire downstream region between the two-dimensional wedge and the shock, and it is supersonic 
since sin (/? - x^j < 1 , as in the classical case. The flow through an attached shock and around a two-dimen- 
sional wedge with ionization, is summarized as follows: 

(1) The shock is a two-dimensional wedge with a shock angle, ft , given by the oblique shock 
curves with ionization for the specified wedge angle, x > and upstream Mach number. 

(2) The streamlines turn by an angle of p - x when they cross the shock, and they flow down- 
stream in parallel with the surface of the wedge. 

(3) The flow variables are supersonic and constant throughout the region between the two- 
dimensional wedge and the shock and along the surface of the wedge. p dn and p dn are given by equa- 
tion (63), while M dn = l/sin^ - xj ■ 

As such, the attached shock and the flow through the shock and along the surface of the two- 
dimensional wedge with ionization are the same as for an attached shock and the flow around a 
two-dimensional wedge without ionization, except for the shock angle and the numerical values of the 
flow variables. 
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5.3.2 Attached Shocks With Ionization Near the Vertex 


Rather than ionizing the entire volume of a shock, one would like to restrict the ionization to a 
volume near the vertex of the two-dimensional wedge to reduce the energy required for ionization. This is 
illustrated in figure 17 where it has been assumed that the ionization level in the volume covering the part 
of the shock near the vertex of the two-dimensional wedge is sufficient to drive the shock intensity to its 
minimum and that there is no ionization around the outer part of the wedge. 



Figure 17. Flow through an attached shock and around a two-dimensional wedge 
with ionization near the vertex. 


The streamlines through both parts of the shock flow in parallel with the surface of the wedge 
without expanding or contracting. Near the surface of the wedge the shock angle, (3 , the angle between 
shock and the streamlines downstream of the shock, f3 - x ■> ar| d the downstream flow variables — 
M dn , /i t [ n , p t [ n , and /7 t | n —are the same as in section 5.3.1, where the ionization covers the entire shock. 
Away from the surface of the wedge the shock angle, [3 , the angle between shock and the streamline 
downstream of the shock, (3- x , and the downstream flow variables— Af t | n , p dn , p dn , and M d n — are the 
same as without ionization. To verify this, recall that the criteria, for the streamlines between the shock and 
the two-dimensional wedge to be compatible with the fluid flow equations, are as follows: 

(1) They must fill the region between the shock and the two-dimensional wedge. 

(2) They must vary continuously downstream of the shock satisfying Bernoullis equation. 

(3) The angle between the shock and the streamline downstream of the shock must be consistent 
with the upstream Mach number and properties of the shock at the point where the streamline crosses the 
shock. 
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The first two items are clearly satisfied by the flow illustrated in figure 17 where the flow variables 
are constant in each streamline since the flow in the streamlines is neither compressive or expansive. Fur- 
thermore, as indicated in the derivation of section 5.3.1, the angle between the shock and the a streamline 
downstream of the shock is only dependent on the upstream Mach number and the properties of the shock 
at the point where the streamline crosses the shock. 

Finally, since the flow variables along the surface of the two-dimensional wedge are the limiting 
values of the flow variables in the streamlines approaching the surface of the wedge, the performance of 
a two-dimensional missile or aircraft modeled by the wedge with ionization near the vertex of the wedge 
will be indistinguishable from the case where the “entire” volume around the shock is ionized. As a mat- 
ter of practice, however, the ionized volume near the vertex should be wide enough to guarantee that the 
entire boundary layer is contained in the region where the ionization effects the flow variables (sec. 6.2). 

5.3.3 Detached Shocks With Ionization 

The flow pattern for the detached case with ionization is shown in figure 18, where an “upper bar” 
is used to denote the variables that change with ionization. Although the shock wave is detached and nor- 
mal to the flow at the vertex of the two-dimensional wedge, it does not stand off from the vertex. Indeed, 
the stand-off distance, r , for a detached shock is 

_ 2 r 

T = -rz =f . (66) 

3|_Pdn/ Pup — lj 


where r is the radius of curvature of the vertex. 32 As such, r = 0 for an ideal wedge. Moreover, the local 
shock angle, [5 (x) , of the detached shock decreases with increasing x, starting at 90° for x = 0 (at the ver- 
tex). As such, the flow along each streamline is different, while the limiting values of the flow variables as 
the streamlines approach the surface of the wedge are required for aerodynamic analysis. 


To analyze the flow immediately behind the shock, in a streamline that crosses the shock at x , 
compute the pressure and density via the modified pressure and density ratios with ionization of equation 
(58), using the Mach number of the upstream flow normal to the shock at x, M” p (x) = M up sin(/I(x)) . 
This yields 

2 2 

PdnM /KpSin^W)) +i p dn (x) (r +1 )( M up sm (/?M)) 

= 7 r and = » . (67) 

Pu P ( 7 + 1 ) Pu P y(M up sin(/3(jc))) +1 
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Figure 18. Flow through a detached shock and around a two-dimensional wedge with ionization. 


As in the attached case, Af dn = 1 with ionization. Unlike the attached case, however, the angle 
between the shock and the streamline immediately behind the shock is unknown. As such, use the fact that 
the shock does not effect tangential velocity; i.e., w dn (v) = p (x) , allowing computation of the tangen- 
tial Mach number immediately behind the shock via 

jr.t ( \ _ M dn ( A ) _ M dn ( A ) ^up _ M up( A ) Qip 
LV1 a„\X) — — — = — = 

Qin ( A ) C up C dn (v) C up C dn (v) 

’ <68) 

where C dn (x) = ^y/J dn (w)/ p dn (x) . Given M dn (x) , one may now compute the angle, 6 , between the 
shock and the streamline crossing the shock at x immediately behind the shock, in terms of the local shock 
angle via 


while 


6 (x) = atan 


*4,M 


= atan 


Qln ( A ) 


C U pM U p cos 


(P( x )) 



( 69 ) 


( 70 ) 
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As such, unlike a classical detached shock, the flow immediately behind the shock in a stream- 
line is supersonic. Furthermore, since the flow in a strea ml ine is compressive downstream of the shock, 
Bernoulli’s equation implies that the Mach number of the flow decreases along the streamline. 16 Since 
the streamline is already downstream of the shock, however, the Mach number of the flow cannot drop 
below 1 (without a second shock). As such, the flow is supersonic over the entire length of the streamline. More- 
over, this analysis is independent of x, and as such, the flow is supersonic over the entire region between the 
two-dimensional wedge and the detached shock. 

Finally, to estimate the flow variables along the surface of the wedge, observe that 


lim A/ dn (x) = lim 
x — >0 x — ^0 



lim , 

x — >0 


1 + 


M u p cos 


{P( x )) 


C, 


\2 


up 


^dn (a) 


= 1 


(71) 


immediately behind the shock since fi(x) = 90° at the x = 0 where the shock wave is locally normal. 
Since the flow cannot drop below Mach 1 anywhere in the streamline, this implies that lim M dn (x) = 1 

x — ^0 

over the entire downstream part of the streamline, while Bernoulli’s equation implies that in the limit as 
x — > 0 , all of the flow variables are constant over the entire downstream part of the streamline. 


Since the shock is locally normal at x = 0 ? one can compute the flow variables in the limiting 
streamlines as x — » 0 , via the modified pressure and density ratios with ionization of equation (60) using 
the actual upstream Mach number: 


, im gdnW = r<+ 1 a„ d 

x ~ >0 Pup (7 + 1) x— >0 p U p 


(r + i)Mu P 

r^up + i 


(72) 


while 


lim M dn (x) = 1 . (73) 

x — ^0 

Finally, since the flow along the surface of the two-dimensional wedge is determined by the limit- 
ing value of the flow in the streamlines as x — > 0 , the flow variables along the surface of the two-dimen- 
sional wedge with ionization are also given by equations (72) and (73). 

The above flow through a detached shock and around a two-dimensional wedge with ionization 
may be summarized as follows: 

(1) The shock is normal at the vertex of the two-dimensional wedge with a stand off 
distance f = 0 . 

(2) The local shock angle, f3(x), is 90° at x = 0 ; i.e., at the vertex of the two-dimensional 
wedge), decreasing with increasing x. 
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(3) The streamlines cross the shock and turn by an angle, 0 (v) , immediately behind the shock, 
continuing to turn downstream until they reach a steady state value. 

(4) The pressure, density and Mach number of the flow in a streamline immediately behind the 
shock are given by equations (57) and (60). 

(5) The Mach number in a streamline decreases downstream of the shock but remains supersonic 
while the pressure and density increase downstream. As such, the flow variables are supersonic through- 
out the region between the two-dimensional wedge and the shock. 

(6) The pressure and density along the surface of the two-dimensional wedge are given by equa- 
tion (72) while the Mach number along the surface of the two-dimensional wedge is 1. 

Note, the above analysis implies that the flow variables along the surface of the two-dimensional 
wedge are the same as those that would result from a two-dimensional analysis with a normal shock 
wave. 
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6. AERODYNAMICS AND PERFORMANCE OF A HYPERSONIC LIFTING BODY 


Over the past decade, the driving force behind much of the research in shock structure analysis 
has been its potential application to aeronautics. 2 ’ 5 ’ 6 ’ 32 ’ 34 ’ 35 Indeed, the reduction in pressure and density 
behind a shock wave propagating through a weakly ionized gas suggest the possibility of reducing both 
wave (pressure) drag and friction drag on the underlying missile or aircraft. Unfortunately, these effects 
may be compromised by the broadening of the shock wave, which increases the Mach number of the nor- 
mal airflow through the shock and the pressure and density behind the shock. 

Indeed, our simulations indicate that ionization may increase drag on a cylindrical missile at zero 
angle of attack. For a lifting body (or a missile/aircraft at a high angle of attack), however, the broadening 
of the shock wave can be beneficial, increasing the lift on the lower surface of the lifting body, without 
requiring an increased angle of attack which would also increase drag. 

The purpose of this section is to use the aerodynamic tools developed in section 5 to analyze the 
aerodynamics and performance of a hypersonic lifting body with ionization at various Mach numbers and 
angles of attack, comparing its performance to a lifting body of the same design without ionization. For 
this purpose a two-dimensional hypersonic lifting body whose geometry is indicated in figure 19 will be 
used. The analysis of the flow along the lower forward surface of the lifting body with ionization is per- 
formed via the techniques developed in section 5. Since the ionization does not extend behind the shock, 
however, one can use classical techniques to analyze the flow along the lower rear and upper surfaces of 
the lifting body. 24 ’ 31 



A summary of an aerodynamic performance analysis of the two-dimensional hypersonic lifting 
body with and without ionization, parameterized by altitude, Mach number, and angle-of-attack appears 
in section 6.1, while the energy required by the ionization system is evaluated, and the trade-offs between 
ionization energy and the aerodynamic performance of the lifting body are reviewed in section 6.2. Finally, 
a summary of the classical two-dimensional aerodynamic analysis techniques employed in the analysis 
appears in appendix C. 
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6.1 Aerodynamic Analysis 


In the following, a Mach 4 lifting body operating at 70,000 ft with the operational and geometric 
parameters indicated in table 3 is considered. The lift over drag (L/D) curves for the lifting body with and 
without ionization at various angles of attack appear in figure 20 (a). Indeed, the L/D curves with ioniza- 
tion are consistently higher than without ionization, peaking at 8.27 at an angle of attack of 2.41° with ion- 
ization, whereas the peak L/D without ionization is 4.21 at a = 5.04° . Note, the “kink” in the L/D curve 
with ionization is at the point where the shock detaches. Since L/D is maximized with an attached shock 
in this example, the “with ionization” curves in figures 20 and 21 are truncated at a = 3.38° , where the 
shock detaches (except for the L/D curve). There are, however, some combinations of Mach number and 
nose cone angle where one may want to consider the possibility of flying with a detached shock, such as 
the Mach 2, 70,000-ft case summarized in table 4. 


Table 3. Parameters of lifting body. 


M up 

Altitude 

V 

(in 

V 

Width 

4 

70,000 

10° 

5 m 

10 m 

1 m 


From the lift and drag curves of figures 20 (b) and (c), it is apparent that the L/D increase with 
ionization is due to increased lift with drag only minimally decreased (table 4). Also note that even though 
the air temperature along the lower forward surface of the lifting body with ionization is twice that with- 
out ionization (fig. 21), the temperature increase in the boundary layer is smaller with ionization since the 
Mach number of the air flow along the boundary layer is smaller, resulting in similar skin temperatures 
with and without ionization as shown in figure 20 (d). Finally, the air flow variables along the lower for- 
ward, lower rear, and upper surfaces of the lifting body are plotted against angle of attack with and without 
ionization in figure 21. 

From the point of view of the aircraft designer, one would prefer to compare the two vehicles 
trimmed to provide the same lift, and hence carry the same payload and fuel fraction, rather than trimmed 
at the point of maximum L/D. As such, assume that the lifting body without ionization is trimmed to maxi- 
mize L/D, as above, while the lifting body with ionization is trimmed to provide the same lift as without 
ionization. To this end, one would trim the lifting body with ionization in the above Mach 4 example, to 
fly at an angle of attack of 0.58°, generating the same 88,300 N of lift as the optimally trimmed vehicle 
without ionization, but with reduced drag and skin temperature (table 4). Although the L/D in this case is 
less than in the optimally trimmed case (7.55 versus 8.27), it allows one to compare the performance of 
alternative aircraft designed to carry the same payload and fuel fraction, with and without ionization. 

Finally to characterize the effect of Mach number and altitude on the performance of the lifting 
body, the key performance parameters for the two-dimensional lifting body of figure 19 with and with- 
out ionization, are tabulated for Mach numbers ranging 2 to 12 and altitudes ranging from 40,000 ft to 
100,000 ft in tables 4 and 5. 
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xIO 5 




(a) (b) 


xIO 4 




(c) (d) 

Figure 20. Performance of a Mach 4 two-dimensional lifting body at 70,000 ft versus 
angle of attack with (solid) and without (dashed) ionization: (a) L/D, (b) lift, 
(c) drag, and (d) skin temperature. 
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xIO 4 











Figure 21. Air flow variables for Mach 4 two-dimensional lifting body at 70,000 ft versus 
angle of attack with (solid) and without (dashed) ionization. 
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Table 4. Performance of a two-dimensional lifting body at 70,000 ft versus Mach number 
with and without ionization. 


Mach 

Parameter 

With 

Ionization 
Max L/D 

With 

Ionization 

=Lift 

Without 
Ionization 
Max L/D 

Parameter 

With 

Ionization 
Max L/D 

With 

Ionization 

=Lift 

Without 
Ionization 
Max L/D 

2 

Angle of attack 

0° 

** 

5.41° 

Shock angle 

Detached 

** 

46° 


L/D 

7.41 

** 

4.61 

Lift 

50,100 N 

** 

47,500 N 


Drag 

6,876 N 

** 

10,297 N 

Skin temperature* 

491 K 

** 

417 K 

3 

Angle of attack 

2.03° 

0.83° 

5.54° 

Shock angle 

55.3° 

44.3° 

32.9° 


L/D 

8.5 

7.93 

4.31 

Lift 

117,900 N 

68,300 N 

68,300 N 


Drag 

13,900 N 

8,610 N 

15,900 N 

Skin temperature* 

686 K 

637 K 

670 K 

4 

Angle of attack 

2.41° 

0.58° 

5.04° 

Shock angle 

44.2° 

36.6° 

27.2° 


L/D 

8.27 

7.55 

4.21 

Lift 

163,000 N 

88,300 N 

88,300 N 


Drag 

19,700 N 

11,700 N 

21,000 N 

Skin temperature* 

1,039 K 

976 K 

1,004 K 

6 

Angle of attack"* 

1.64° 

0° 

4.14° 

Shock angle 

35.3° 

30° 

21.9° 


L/D 

8.2° 

7.34 

4.19 

Lift*** 

241,900 N 

141,000 N 

134,000 N 


Drag 

29,500 N 

19,200 N 

32,300 N 

Skin temperature* 

2,030 K 

1,941 K 

1,925 K 

8 

Angle of attack 

1.06° 

** 

3.4° 

Shock angle 

31.5° 

** 

19.2° 


L/D 

8.31 

** 

4.16 

Lift 

341,800 N 

** 

182,700 N 


Drag 

41,100 N 

** 

44,000 N 

Skin Temperature* 

3,404 K 

** 

3,192 K 

10 

Angle of attack 

1.02° 

** 

2.75° 

Shock angle 

30.5° 

** 

17.5° 


L/D 

8.44 

** 

4.2 

Lift 

513,500 N 

** 

235,600 N 


Drag 

60,800 N 

** 

56,100 N 

Skin temperature* 

5,210 K 

** 

4,812 K 

12 

Angle of attack 

0.6° 

** 

2.11° 

Shock angle 

28.7° 

** 

16.2° 


L/D 

8.56 

** 

4.26 

Lift 

645,000 N 

** 

287,300 N 


Drag 

75,400 N 

** 

67,400 N 

Skin temperature* 

7,334 K 

** 

6,788 K 


* Measured at the lower forward surface of the lifting body. 

'* Requires a negative angle of attack, which is not supported by our code. 
'* Approximation for the equal lift case. 


Table 5. Performance of a Mach 4 two-dimensional lifting body versus altitude with 
and without ionization. 


Alt (ft) 

Parameter 

With 

Ionization 
Max L/D 

With 

Ionization 

=Lift 

Without 
Ionization 
Max L/D 

Parameter 

With 

Ionization 
Max L/D 

With 

Ionization 

=Lift 

Without 
Ionization 
Max L/D 

40,000 

Angle of attack* 

2.22° 

0.26° 

4.94° 

Shock angle 

43.2° 

35.5° 

27.1° 


L/D 

8.44 

7.41 

4.30 

Lift 

600,800 N 

304,100 N 

304,100 N 


Drag 

71,200 N 

41,000 N 

79,000 N 

Skin temperature* 

963 K 

905 K 

939 K 

70,000 

Angle of attack* 

2.41° 

0.58° 

5.04° 

Shock angle 

44.2° 

36.6° 

27.2° 


L/D 

8.27 

7.55 

4.21 

Lift 

163,000 N 

88,300 N 

88,300 N 


Drag 

19,700 N 

11,700 N 

21,000 N 

Skin temperature* 

1,039 K 

976 K 

1,004 K 

100,000 

Angle of attack* 

2.65° 

0.72° 

5.33° 

Shock angle 

45.6° 

36.9° 

27.5° 


L/D 

8.06 

7.3 

4.1 

Lift 

44,700 N 

23,500 N 

23,500 N 


Drag 

5,530 N 

3,220 N 

5,720 N 

Skin temperature* 

1,094 K 

1,021 K 

1,043 K 


* Measured at the lower forward surface of the lifting body. 
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Note that a negative angle of attack would be required to match the lift with ionization to the 
optimal lift without ionization in the Mach 2, 8, 10, and 12 cases at 70,000 ft. Since this would not take 
advantage of the available lift with ionization and is not supported by our code, an “equal lift” analysis 
was not performed in these cases, indicated by “**” in table 4. Also note that some type of active cooling 
system would be required to deal with the skin temperatures for the Mach 8-12 cases. 

6.2 Energy Requirements 

Since the energy used by the ionization system could (theoretically) have been used to increase the 
engine thrust, it should be incorporated into the L/D analysis to properly evaluate the performance of the 
ionization system on the lifting body of section 6.1. To this end, the “effective” drag of the energy used by 
the ionization system has been included in the L/D and drag of a lifting body with ionization (trimmed to 
provide the same lift as an optimally trimmed lifting body without ionization), and it has been compared 
to the L/D and drag without ionization. 

To estimate the energy requirements for an ionization system, assume that ionization is only used 
near the vertex of the two-dimensional nose cone of the lifting body (fig. 22). As indicated in section 5.3.2, 
the shock in the ionized region takes the angle determined by the oblique shock curves with ionization 
of section 5.2, and the flow variables in the streamlines through the ionized part of the shock take on the 
values determined by the modified pressure and density ratios with ionization of section 5.2. Similarly, the 
shock angle in the region without ionization and the flow variables in the streamlines passing through that 
part of the shock are characterized by the classical oblique shock curves and pressure and density ratios. 
Furthermore, both sets of streamlines flow in parallel with the surface of the lifting body, while the flow 
along the surface of the lifting body is determined by the streamlines that “pass through” the boundary 
layer. As such, it would be theoretically possible to limit the ionized region to the thickness of the bound- 
ary layer, 1 mm or less at the far end of the lower rear surface of the lifting body. In practice, however, it 
is unlikely that one could focus an ionization source that precisely. As such, for the purpose of the pres- 
ent analysis, assume that the ionized part of the shock is 1.5 cm across (a half wavelength for an X-Band 
microwave source). Since the lifting body is 1-m wide, this implies that the source is required to ionize an 
area of 1.5x10“- m 2 . 



Figure 22. Shock, flow, and boundary layer with ionization near the nose of the lifting body. 
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The minimum power required to drive the shock intensity to its minimum value for the five flight 
regimes with equal lift with and without ionization (tables 4 and 5), was evaluated using the shock struc- 
ture analysis algorithm of section 4, and is summarized in table 6 together with the equivalent drag force 
(power/upstream velocity). The effective drag with ionization is then taken to be the sum of the actual 
drag with ionization plus the equivalent drag of the power required by the ionization source. The effective 
L/D and drag reduction are then computed, both of which are reduced by the equivalent drag of the power 
required by the ionization source. Finally, following Beaulieu et al., define the energy effectiveness of the 
ionization system as the ratio of the effective drag reduction divided by the equivalent drag of the power 
required by the ionization source. 34 These data are summarized in table 6 for the five flight regimes with 
equal lift with and without ionization. 


Table 6. Power required for ionization and effective drag with ionization 
for the lifting body at 70,000 ft. 


Flight Regime 

Mach 3 
(70,000 ft) 

Mach 4 
(70,000 ft) 

Mach 6 
(70,000 ft) 

Mach 4 
(40,000 ft) 

Mach 4 
(100,000 ft) 

Power to ionize 1 .5x1 0 -2 m 2 of the shock 

2.42x10 3 W 

5.22x10 3 W 

1.08x10 4 W 

1.94x10 4 W 

2.45x10 s W 

Equivalent drag force 
(Power/ upstream velocity) 

2.72 N 

4.4 N 

6.208 N 

16.9 N 

202 N 

Drag with ionization 
(table 4) 

8,610 N 

11,700 N 

19,200 N 

41,000 N 

3,220 N 

Effective drag with ionization 
(drag with ionization + equivalent drag) 

8,612 N 

11,704 N 

19,206 N 

41,017 N 

3,422 N 

Effective L/D 

7.931 

7.544 

7.341 

7.414 

6.86 

Effective drag reduction 
(drag without ionization -effective drag 
with ionization) 

7,288 N 

9,296 N 

13,094 N 

37,983 N 

9,142 N 

Energy effectiveness 
(effective drag reduction /equivalent drag 
of ionization power) 

2.68x10 3 

2.11 xIO 3 

2.11 xIO 3 

2.25 xIO 3 

4.53X10 1 


With the exception of the Mach 4, 100,000-ft case, the energy effectiveness of the ionization sys- 
tem exceeds 2,000, indicating that an ionization system would significantly increase the performance of 
the lifting body, even if one were to increase the power required by the ionization system by an order of 
magnitude (or two) to incorporate thermal and volumetric inefficiencies in the analysis. Indeed, such a 
system would achieve a net performance gain even in the 100,000-ft case, though with reduced energy 
effectiveness. Finally, note the nonmonotonicity of the data with altitude, which is consistent with the 
observations based on the simulations of section 4. 
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7. CONCLUSIONS 


In this section the assumptions underlying the electrofluid dynamics model are reviewed, a qualita- 
tive comparison of the above analysis and simulations with the available experimental data is provided, 
and the potential applications of the shock structure modification technology in aeronautics are com- 
mented upon. 


7.1 Electrofluid Dynamics Model 

The driver for the shock structure analysis is the electrofluid dynamics model of equations (1 1)— (20). 
Indeed, this model is composed of 6 classical conservation laws (for either the entire weakly ionized gas 
or the individual species), plus Gauss’ Law, and the required auxiliary equations. The remaining conserva- 
tion laws are implied by the assumptions that u l = u n , T l = T n , and that the electron momentum is negli- 
gible. As such, the efficacy of the resultant analysis is dependent on these assumptions and the remaining 
assumptions summarized in table 2. 

Of equal significance are the assumptions that are not made, as follows: 

• Not assuming quasi-neutrality and including the resultant electric field and electrostatic body forces/ohmic 
heating in the model. 

• Treating the electron velocity as an independent variable. 

• Including viscosity in the model. 

Although, to our knowledge, all previous shock structure analysis models have assumed that 
u e = u n = u , the electron velocity may increase by several orders of magnitude in the electric field across 
the shock, significantly impacting the shock structure. 1 ’ 9 ’ 12 ’ 13 In particular, the increased electron veloc- 
ity decreases the electron drift length while increasing the electron body force/ohmic heating term in the 
energy equation, making it the primary coupling mechanism between the ionization and the fluid flow. 

Similarly, the viscosity in the fluid dynamics model significantly influences the shock structure. 
In particular, the jump in electron temperature across the shock is driven by the pressure jump, which is 
independent of shock width, but the loss mechanisms, which reduce the electron temperature jump in the 
shock, are dependent on the shock thickness. Furthermore, the viscosity of the flow dampens the rapid 
decrease of the flow intensity as it converges to its minimum value downstream of the shock and enhances 
the decrease of the Mach number upstream of the shock. 

Finally, note that the thermal radiation term, Q, can be used to add a thermodynamic model to the 
electrofluid flow model. Such a model might include thermal radiation, if significant, an external heat 
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source, or a feedback system in which Q is dependent on the fluid flow variables. Since Q is additive to 
Q e , it will either add or subtract from the ohmic heating produced by the electrodynamics model, but 
should not fundamentally modify the shock structure. 

7.2 Comparison With Experiment 

The experimental literature on ionized gas/shock wave interaction can be divided into the 
following: 

• Macroscopic effects— where one observes the effect of the ionization indirectly, via its effect on the 
shock wave. 

• Microscopic effects— where one observes the effect of the ionization directly via its effect on the fluid 
flow and electrodynamic variables across the shock 

Unfortunately, the available electron temperature and density measurements, are not sufficiently 
precise to permit a quantitative comparison with the simulations of section 4, where a 0.01-eV change in 
electron temperature can fundamentally effect the shock structure. As such, this review is restricted to a 
qualitative comparison. 

The most common macroscopic effects include the following: 3-8 ’ 32 ’ 36 

• Reduction in shock intensity. 

• Broadening of the shock wave. 

• Increased standoff distance for a detached shock wave. 

• An enhancement of the observed effect when an external electric field is applied in the direction oppo- 
site to the flow, and a degradation of the observed effect when an external electric field is applied in the 
direction of the flow. 5 ’ 37 

• The apparent complete dispersion of the shock wave with sufficiently high ionization levels. 5 ’ 7 ’ 14 

• Asymetric time constants associated with the change in shock structure (changing in microseconds 
when the ionization is turned on, but requiring milliseconds to revert back to the original shock structure 
when the ionization is turned off). 38 

The first three items are due to the reduced density ratio across the shock in the presence of ion- 
ization. Indeed, these observations are all made with either Schlieren or shadowgraph images, which are 
proportional to either the derivative of the second derivative of the density, 39 and as such, the intensity of 
these images is reduced when the density ratio across the shock decreases in the presence of ionization. 

The observed broadening of the shock wave is also attributable to the reduced density ratio across 
the shock in the presence of ionization via the equality 
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1 


( 74 ) 


Pdn _ 

Pup cot (/J) tan (/J - x) ’ 


which is used in the derivation of the oblique shock curves, 31 where fi is the shock angle and % is the 
nose cone angle. Similarly, the shock standoff distance increases with decreasing density ratio (eq. (66)). 

Since the electric field that drives the changes in shock structure in the above analysis (via the 
ohmic heating term) is always negative; i.e., in the direction opposite to the flow, an externally imposed 
electric field in the direction opposite the flow will enhance its effect, whereas an externally imposed elec- 
tric field in the direction of the flow will detract from the effect of the ionization. 

Although the apparent complete dispersion of the shock wave in the presence of ionization would 
appear to contradict the fact that the shock intensity cannot be driven below a minimum value; 3 ' 7 ' 14 
the Schlieren and shadowgraph diagnostics used in these experiments image the derivative and second 
derivative of the gas density, respectively. 39 As such, they are dependent on the distance upstream and 
downstream of the shock over which the ionization interacts with the fluid flow as well as the density 
ratio. Indeed, as the ionization level increases in these simulations, the distance upstream and downstream 
of the shock over which the ionization interacts with the fluid flow continues to increase after the shock 
intensity and flow density have reached minimum values. As such, the (average) amplitude of the Schlie- 
ren and shadowgraph images will continue to decrease with increasing ionization levels. This is illus- 
trated in figure 23 for a Mach 2 shock at 70,000 ft, which roughly approximates the environment where 
the experimental observations were made. Although the downstream density increases slightly for the 
1.4-eV case (probably numerical error for this example, which is on the fringes of instability), the distance 
over which the density jumps across the shock is about 40 pm for the 1-eV case and about 150 pm for the 
1.4-eV case. As such, depending on the spacial resolution of the diagnostic systems employed the (aver- 
age) amplitude of the Schlieren and shadowgraph images for the 1.4-eV case will be about a third of that 
for the 1-eV case. 




Figure 23. Density for a Mach 2 shock at 70,000 ft with (a) 1-eV and (b) 1.4-eV sources. 
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Finally, note that the asymmetric time constants observed experimentally, 38 cannot be modeled by 
our steady state electrofluid dynamics model. Given the complex feedback loops coupling the fluid flow 
and electrodynamic variables in the model of figure 4, however, such an effect would not be unreason- 
able. 


Unlike the macroscopic effects, high (spatial) resolution instrumentation is required to directly 
observe the microscopic effects of ionization on the shock structure. To increase the spatial scale of the 
region around the shock wave, most researchers work at low pressures (typically 3-5 Torr), thereby facili- 
tating the direct observation of the fluid flow and electrodynamic variables within the shock wave. Liu, 
Whitten, and Glass; and Liu and Glass use a Mach-Zehnder interferometer to observe the structure of Mach 
13-16 shock waves in 5 Torr argon, in a 10-cm by 18-cm. hypervelocity shock tube. 12 ’ 40 In these experi- 
ments they observe a “small” electron temperature peak at the shock and a factor of 2 jump in electron 
density behind the shock wave. In a similar series of experiments in Mach 17-21 neon, McIntyre et al., 
saw an order-of-magnitude or more jump in electron density at the shock and observed the increased ion- 
ization rate at the shock indirectly via measurements of the excited state populations. 10 Moreover, they 
saw an overshoot in the heavy particle concentration behind the shock, which is consistent with the over- 
shoot in our simulations. Finally, Chan and Hershkowitz saw an electron temperature peak at the shock 
and the associated electronic double layer. 41 

In a recent paper, Bletzinger, Ganguly, and Garscadden reported on an experimental program in 
3-5 Torr nitrogen, performed in a 3-cm-diameter shock tube with an electric shock driven Mach 1.76 
shock wave propagating through a 20 ma weakly ionized gas. 35 They summarized the results of these 
experiments as follows: 

[ W \ e have provided direct evidence for the large jump in ionization potential, the fast 
increase of electron excitation temperature as indicated by the sudden increase in emission rate, 
and the increase in the emission rate and the increase in n e caused by the formation of strong 
double layers, at the shock front, in non-equilibrium ionizations. One of the important implica- 
tions of the strong- double-layer-induced local electron heating is that a very localized energy 
dispersion and gas heating at the shock front is possible so that shock wave dispersions and 
modifications of the shock jump conditions cannot always be determined from the steady state 
or spatially averaged gas properties . 35 

In particular, Bletzinger et al., observed the following: 35 

• A negative electric field across the shock. 

• An increased ionization rate at the shock. 

• The formation of a double layer at the shock. 

• A drop in the density ratio across the shock. 

Based on these observations, Bletzinger et al., concluded that the electron temperature peaks at the 
shock and that the observed shock wave dispersion is due to local ohmic heating at the shock, induced by the 
double layer. As such, the microscopic shock structure model resulting from the Bletzinger et al. experiments 
is qualitatively identical to the model predicted by the above described shock structure analysis and 
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simulations. 35 The experimental and theoretical models differ quantitatively in certain respects — a 10-eV 
electron temperature jump in the experimental model versus <1 eV in the analytical model and a 
150 percent electron density ratio in the experimental model versus more than an order of magnitude in 
the analytical model. 


7.3 Application in Aeronautics 

The potential applications of the shock structure modification in a weakly ionized gas are depen- 
dent on the trade-offs between its various effects, as follows: 

• Reducing the pressure and density jump across the shock. 

• Broadening of the shock wave. 

• Increasing the temperature jump across the shock. 

For the lifting body of section 6, the broadening of the shock wave is the primary source of the 
increased L/D with ionization— increasing the Mach number of the normal air flow through the shock 
and hence the pressure below the lifting body at a lower angle of attack than without ionization, thereby 
increasing lift without a commensurate increase in drag. Moreover, in this application, the skin tempera- 
ture of the lifting body is essentially unchanged, since the decreased downstream Mach number with ion- 
ization compensates for the increased downstream air flow temperature (eq. (97)). 

The above analysis would also apply to a missile or aircraft flying at a high angle of attack, or the 
flow in a duct or inlet, where a broad shock angle is also beneficial. Moreover, for a “blunt” nosed missile 
or aircraft, where the shock angle is already “broad,” the additional broadening of the shock wave with 
ionization will only minimally increase the Mach number of the normal air flow while decreasing the pres- 
sure and density behind the shock with a commensurate decrease in wave and friction drag. For a “fine” 
nosed missile or aircraft, however, the broadening of the shock wave with ionization will counteract the 
reduction of the pressure and density jumps across the shock, possibly increasing drag. 
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APPENDIX A- ADDITIONAL SHOCK STRUCTURE ANALYSIS DATA 


The purpose of this appendix is to expand on section 4.1 by providing detailed shock structure 
analysis plots for some additional cases (figs. 24-26). 
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A.l Shock Structure of a Mach 4 Flow at 40,000 ft With a 1.4-eV Source 
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Figure 24. Shock structure for a Mach 4 flow at 40,000 ft with a 1.4-eV source: (a) Pressure, 

(b) density, (c) velocity, (d) Mach number, (e) speed of sound, (f) temperature, 

(g) electric field, (h) electron velocity, (i) ohmic heating, (j) ion and electron densities, 
(k) ionization rate, (1) electron temperature. 
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A.2 Shock Structure of a Mach 12 Flow at 70,000 ft With a 0.85-eV Source 
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Figure 25. Shock structure for a Mach 12 flow at 70,000 ft with a 0.85-eV source: (a) Pressure, 
(b) density, (c) velocity, (d) Mach number, (e) speed of sound, (f) temperature, 

(g) electric field, (h) electron velocity, (i) ohmic heating, (j) ion and electron densities, 
(k) ionization rate, (1) electron temperature. 
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A.3 Shock Structure of a Mach 2 Flow at 10,000 ft With a 0.9-eV Source 
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Figure 26. Shock structure for a Mach 2 flow at 10,000 ft with a 0.9-eV source: (a) Pressure, 

(b) density, (c) velocity, (d) Mach number, (e) speed of sound, (f) temperature, 

(g) electric field, (h) electron velocity, (i) ohmic heating, (j) ion and electron densities, 
(k) ionization rate, (1) electron temperature. 
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APPENDIX B- APPROXIMATE ANALYTIC DERIVATION OF THE MODIFIED PRESSURE 

AND DENSITY RATIOS 


The purpose of this appendix is to provide an alterative approximate analytic derivation of the 
effective Mach numbers that bypasses the empirical derivation used in section 5.1. To this end, start with 
the inequalities 


( y + 1 X and 

Pup ~ (y-l)M„ p + 2 ~ Pa p- (r + 1) ' 


(75) 


where Pdn/Pup and Pdn/Pup are u PP er bounded by the “no ionization” case since ionization reduces 
the density and pressure ratios across the shock, and lower bounded by the “no shock” case where 
Pdn = p U p and p t [ii = Pup ■ Substituting the density inequalities into the expression for the downstream 
shock intensity 


-\j{ PdnPdn )/( Pup Pup ) — ^up 

with ionization, one may derive tighter upper and lower bounds for p^ n / p up via 

m 2 = PdnPdn < Pdn {7 +l ) M up 
Pup Pup Pup (y-l)M^p+2 

or equivalently 

Pdn >M 2 (r-l)M u 2 p + 2 _ (y-l)M u 2 p+ 2 

Pup” up (r + i)M 2 p - (r+i) ‘ 


(76) 


(77) 


(78) 


Similarly, one may derive an upper bound for Pdn/Pup using the opposite inequality in equation (75) 
yielding 

m 2 = PdnPdn > Pdn j = Pdn (79) 

Pup Pup Pup Pup 

Although there remains a considerable gap between the upper and lower bounds of equations (78) 
and (79), they are tighter than the bounds of equation (75). In particular, both the upper and lower bounds 
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grow with Myp for large upstream Mach numbers. As such, approximate Pdn/Pup by a value that is 
midway between the upper and lower bounds of equations (78) and (79); i.e., 


2 (y-l)Mj? p + 2 

Pdn up (r + 1) r<, + ! 

Pup ^ (y + l) 


(80) 


Although equation (80) has a maximum percent error of 100 /y (74 percent in air and 60 percent 
in a monotonic gas), it represents the “best” estimate for p dn / Pup given the available data. Substituting 
equation (80) into equation (76) yields a corresponding approximation for the density ratio 

Pdn Jr + iKp (81) 

Pup yA/ up +l 


Finally, one may obtain expressions for the pressure effective and density effective Mach num- 
bers by equating the expressions for Pdn/Pup an d p dn / p up of equations (80) and (81) with the defining 
equalities for M up and M up of equation (53): 


yielding 


while 


2 r^u P - (r - ! ) = Pdn ^ 7^u P + 1 and (r + i)^up _p dn „(r + 1)^ 


(y + !) 


1 up 


Pup (7 + 1 ) (y - 1) A7 2 + 2 Pup yM 2 + 1 


(82) 
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and M 2 


1 + M 
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(83) 
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f 2 M 2 ^ 

Z.iK/ 0 p 


1+M 


tK P = ■ 


up ; 


(84) 


As such, with the aid of the (admittedly “loose”) approximation of equation (80), one can replicate 
the theory of section 5.1 — including the empirical formula for the effective Mach number with ionization. 
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APPENDIX C -TWO-DIMENSIONAL AERODYNAMIC ANALYSIS WITHOUT IONIZATION 


The purpose of this appendix is to provide an overview of the classical two-dimensional aerody- 
namic analysis techniques used in section 6 including the following: 

• The characterization of the shock and the flow behind the shock without ionization, 

• The Prandtl-Meyers expansion analysis used to evaluate the flow on the top and lower rear surfaces of 
the two-dimensional lifting body. 

• The wave drag and lift formulae. 

• The model used for evaluating the friction drag due to the turbulent flow along the surfaces of the lift- 
ing vehicle. 

• A method for computing skin temperature of the lifting body. 

• The boundary layer thickness. 


C.l Shock Angle and Downstream Flow 

For a hypersonic lifting body without ionization, one may assume that the shock is attached— in 
which case, its properties are similar to the properties of the attached shock with ionization described in 
section 5.3, using the classical pressure and density ratio formulae and oblique shock curves rather than 
the modified formula of section 5.3. In particular, the attached shock is wedge shaped, the streamlines 
downstream of the shock are parallel to the surface of the wedge, and the flow variables are constant in the 
region between the shock and the wedge. As such, one may 

(1) Compute the shock angle using the classical oblique shock curves of figure 15 or by solving 


M 


up 


-J-l 2 

sin (J3) \j (y+ 1) cot (P) tan (p-%)- (y- 1) 


(85) 


for p given the wedge angle, / , and M up 

(2) Use the classical pressure and density ratios of equation (52) to compute the pressure and 
density downstream of the shock. 
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via 


(3) Use the inverse of the density ratio to evaluate the normal velocity of the downstream flow 


_P up 

up — 

Pdn Pdn 


n P up n 

M dn = —U 


Sin 


(P)' 


up 


( 86 ) 


(4) Let u dn = u dn /tm((3 - %) ■ 


C.2 Prandtl-Meyers Expansion 


To evaluate the flow on the upper and lower rear surfaces of the lifting body, use a classical 
Prandtl-Meyers expansion since there is no ionization involved in either case. 31 To compute the flow vari- 
ables, M t p, p t p, p tp , and w tp , along the upper surface of the lifting body given the upstream flow variables 
M up, p u p , p up , and i/ up , use the following procedure: 31 

(1) Compute the Mach number of the flow along the upper surface of the lifting body via 


where 


M t p v [oc + v^M U p 



(87) 


( 88 ) 


is the Prandtl-Meyers function plotted in figure 27. 



Figure 27. Prandtl-Meyers function. 
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(2) Given M tp , one may now compute p tp , p tp , and w tp via the total pressure and density 
expressions 


Ptp _ 
Pup 


2 + (y-i)(M U p)' 
2 + (r-l)(M tp ) 2 


r/(r-i) 


( 89 ) 


while 


Ptp 

Pup 


1 

to 

+ 

1 

(M up ) 2 

2 +(r- 1) 

Kf 


-iV(r-i) 


(90) 


M tp ‘^tp tJyp tp / Ptp ’ (91) 

yielding a complete description of the air flow variables along the upper surface of the lifting body. 

Similarly, one can compute the flow variables, A/ lr , p lr , p lr , and iq r along the lower rear surface 
of the lifting body given the flow variables downstream of the shock along the lower forward surface 
of the lifting body, M dn , p dn , p dn , and w dn , via a Prandtl-Meters expansion. Since the ionization at the 
shock does not effect the expansion process, one can use the same Prandtl-Meyers algorithm with or with- 
out ionization. The ionization at the shock, however, indirectly effects the expansion via its effect on the 
flow variables downstream of the shock along the lower forward surface of the lifting body. 

C.3 Wave Drag and Lift 

The wave drag and lift forces on the lifting body are the horizontal and vertical components of 
the forces on the lifting body due to the air pressure on its surfaces, respectively, as shown in figure 28. 
The wave drag and lift forces are computed by multiplying the pressure on the surface by the area of the 
surface projected in the horizontal or vertical direction, respectively. For instance, the lift on a unit width 
of the lower forward surface is L dn = p dn lj cos (a + /]) . Of course, the downward force due to the pres- 
sure on the upper surface of the lifting body, L tp , cancels part of the lift from the pressure on the lower 
surfaces, while the sum of the lift terms must equal the gravitational forces on the lifting body to maintain 
a constant altitude. Similarly, the pressure on the upper and lower rear surfaces of the lifting body partially 
compensates for the wave drag on the lower forward surface, while the sum of the wave drag terms must 
be balanced by thrust from the engine. Finally, to maintain trim, the moments produced by all of the forces 
on the lifting body must be balanced by the moments produced by the control surfaces or attitude control 
jets of the lifting body (not shown in fig. 28). 
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Figure 28. Wave, lift, and friction drag forces on a two-dimensional lifting body. 


C.4 Friction Drag 

To model the friction drag along the surfaces of the lifting body at hypersonic speeds, one must 
account for the turbulent flow along the surface and the effects of the boundary layer. Since the analysis 
requires consideration of the air flow within the boundary layer along the surface of the lifting body, use 
the superscripts “w”, and “m” to distinguish between the air flow variables at the wall or skin of the 
lifting body, the outer edge of the boundary layer, and the mean value between the wall and the outer 
edge, where the air flow variables at the outer edge of the boundary layer coincide with the free air flow 
variables generated by the shock and expansion analysis; i.e., p^ n = p dn , etc. In the following analysis, 
the subscript is dropped when dealing with a generic location. 

Typically a supersonic flow over a flat surface is laminar at the upstream end of the surface, transi- 
tioning to turbulent flow as it passes over the surface, with the friction due to the turbulent flow typically 
dominate. As such, to simplify the friction drag analysis, assume that the flow over the entire surface of 
the lifting body is turbulent. Following Schlichting, the force on a surface of length, /, and unit width due 
to friction drag (in the direction of the flow in parallel with the surface) is given by (ref. 42) 

F = 3.6x10 ~ 6 p m l( M “lp m /p W )”' /5 , (92) 

where the selection of a mean reference level in the boundary layer is based on experimental evidence. 42 
Here, p is the air viscosity, which may be characterized by Sutherland’s equation: 23 

n = 1.458 Xl0“ 6 r 3 / 2 (r + 110.4) _1 . (93) 

Since the air flow variables at the edge of the boundary layer can be evaluated via the shock and 
expansion analysis described above, p m and p m , in equation (90), may be expressed in terms of the 
variables at the edge of the boundary layer. To this end, let 
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where the last equality in equation (92) follows from the observation that the air pressure is constant in the 
boundary layer together with the perfect gas law. 42 Now, following Schlichting 42 



(95) 


where M°° is the Mach number of the flow at the edge of the boundary layer (equal to the Mach number 
of the free air flow). As such, 



(96) 


as required. Similarly, to express /r m in terms of jA°° , observe from Sutherland’s equation (eq. (91)) that, 
at the relatively high temperatures in the boundary layer, ji can be approximated by 1.458xl0 -6 r 1/2 (where 
any approximation error will be further damped out by the fact that viscosity appears to a 1/5 power in the 
friction drag equation (eq. (90)). As such, 
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(97) 


Substituting equations (96) and (97) into the friction drag equation (eq. (92)) now yields the desired 
expression for force due to skin friction on a surface of length, 1, and unit width: 
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(98) 


where <2°° = p°° ^w°° j 2 is the dynamic pressure of the air flow at the edge of the boundary layer and 
R°° = u°° l p°° I p°° is the Reynolds number of the air flow at the edge of the boundary layer. Alternatively, 
upon grouping variables and expressing the viscosity via Sutherland’s equation, F reduces to 


t/ \4/5 / \9/5 Air/ \ 3/10 

F = 3.98 x 10 _3 (p°° ) I u °° ) l 4 / 5 7°°' 


(r°° + 110.4) 1/5 2 + — (m“) 


7/10 


(99) 


Using this model, with the superscript “ 00 ” variables equated to the appropriate free air flow vari- 
ables, one can evaluate the friction drag on each of the three surfaces of the two-dimensional hypersonic 
lifting body of figure 27. 


C.5 Skin Temperature 

To evaluate the thermal loads on the two-dimensional lifting body, one needs to compute the skin 
or wall temperatures, T w . To this end, employ the relationship (ref. 42). 
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C.6 Boundary Layer Thickness 


Following Landau and Lifshitz, 16 the thickness of a boundary layer along a surface at a distance, 
x, from the upstream edge of the surface is 8 = fix /u°° . Taking /A to be the mean value of the viscosity 
in the boundary layer, and expressing fi m via equation (97), this reduces to 


1 



( 101 ) 


Finally, using the “high” temperature approximation of Sutherland’s equation, 

f\ 1 /O 

jl° = 1.458 x 1 0 T , one may express the boundary layer thickness at x via 


8 = 1.435 x 10 -3 


1 



( 102 ) 


In particular, the maximum thickness of the boundary along a surface of length, /, is 


^max 


1.435 xl0“ 3 




(103) 
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